Crop a polygon on svs image

[snigdhaAgarwal]

Is there a possibility of using pyvips to crop a high resolution polygon from an svs image. Something like this but without any pieces of surrounding tissue:

I would like to retain the high resolution of the image by saving the output of the crop as a tiff or svs.

[jcupitt]

Hi @snigdhaAgarwal, do you mean crop a rectangle? Sure, just:

image = pyvips.Image.new_from_file("huge.svs")
area = image.crop(10000, 12000, 5000, 5000)
area.write_to_file("tiny.tif")

Where the four numbers are left / top / width / height in pixels.

[snigdhaAgarwal]

Sorry for the confusion. I meant to say that I have been able to achieve a rectangle but I instead want a polygonal crop according to the shape of the tissue. Right now the rectangle includes pieces of other tissues which I don't want.
I also don't want to lose any resolution. I noticed that opening using new_from_file only loads 1 layer out of the 4 layers of the svs image.

The 4 layers of the same svs when opened with Openslide:

store = OpenSlideStore(path)
grp = zarr.open(store, mode="r")
multiscales = grp.attrs["multiscales"][0]
pyramid = [
    da.from_zarr(store, component=d["path"]) for d in multiscales["datasets"]
]

[jcupitt]

You can use svgload to draw a mask from a set of datapoints, then use the mask to set the image alpha.

Unfortunately, svgload is limited to 32,000 x 32,000 pixels, so it might not go large enough, depending on the size of the area you need to extract. You might need to crop and join in sections.

Yes, libvips always just loads the highest resolution layer of a pyramid. The lower layers are recomputed on save.

[snigdhaAgarwal]

If the other layers are recomputed then why is there a loss in resolution. The original image:

The cropped image:

Notice the pixelation at more zoomed in levels. Is it because I'm replacing the alpha channel(4th band) with the svgload mask?

[jcupitt]

I would guess there's some problem with the viewer. You'd need to give some details about the viewer, processing, parameters, source file, etc.

[snigdhaAgarwal]

Original data file: https://drive.google.com/file/d/16qX9ZLUXm1Dvfd1YCQF-fJNLo3rpSKbc/view

This is the code to get the crop:

y = (6329, 18108,38545,41951,17541,5903)
x = (80269,69341,60116,73031,87933,84101)
points = ([80269, 6329],
 [69341, 18108],
 [60116, 38545],
 [73031, 41951],
 [87933, 17541],
 [84101, 5903])
left = min(x)
top = min(y)
image = pyvips.Image.new_from_file('/Users/snigdha.agarwal/Downloads/5326-1.svs')
crop = image.crop(left, top, max(x) - left, max(y) - top)
svg = f"""
        <svg viewBox="0 0 {crop.width} {crop.height}">
            <polygon style="fill: white; stroke: none" points="
    """
svg += " ".join([f"{x - left}, {y - top} " for x, y in points])
svg += """
            "/>
        </svg>
    """
im=pyvips.Image.svgload_buffer(bytes(svg, "ascii"))[3]
d = crop[:3]
dim = d.bandjoin(im)
dim.write_to_file('trial2.tiff')

I wonder if its a viewer issue because when I try to open this image and see the number of frames, it just shows 1 frame - not a multi-page tiff.

from PIL import Image
img = Image.open('trial2.tiff')

I viewed svs file using the Openslide code from previous comment in a custom version of the napari app. The tiff file can be viwed using graphicConverter 11.

[jcupitt]

libvips saves as a simple flat TIFF by default. Perhaps your viewer is doing nearest neighbour for zoom out?

Try:

dim.write_to_file("x.tif", pyramid=True, tile=True, compression="jpeg")

Deepzoom works well too, if you have a viewer for that.

vipsdisp is a gtk3 viewer that works for most of these formats, though it'll also do nearest neighbour if there's no image pyramid.

https://github.com/jcupitt/vipsdisp

[snigdhaAgarwal]

This worked thanks!

[snigdhaAgarwal]

An aside from the main topic, but is there a way to make this command any faster?

libvips saves as a simple flat TIFF by default. Perhaps your viewer is doing nearest neighbour for zoom out?

Try:

dim.write_to_file("x.tif", pyramid=True, tile=True, compression="jpeg")

Deepzoom works well too, if you have a viewer for that.

vipsdisp is a gtk3 viewer that works for most of these formats, though it'll also do nearest neighbour if there's no image pyramid.

https://github.com/jcupitt/vipsdisp

[jcupitt]

Sorry, that's probably max speed.

I guess you're using libjpeg-turbo?

[snigdhaAgarwal]

I installed libjpeg-turbo but not sure how to ensure that pyvips uses it for the write_to_file command. Right now 1 write takes 2 minutes.

[jcupitt]

That's probably all you need to do to enable it.

You could try changing these lines in svgload.c:

        if( vips_foreign_load_svg_parse( svg, t[0] ) ||
                vips_image_generate( t[0], NULL,
                        vips_foreign_load_svg_generate, NULL, svg, NULL ) ||
                vips_tilecache( t[0], &t[1],
                        "tile_width", VIPS_MIN( t[0]->Xsize, RSVG_MAX_WIDTH ),
                        "tile_height", 2000,
                        "max_tiles", 1,
                        NULL ) ||
                vips_sequential( t[1], &t[2],
                        "tile_height", 2000,
                        NULL ) ||
                vips_image_write( t[2], load->real ) )
                return( -1 );

Change the 2000 to something like 10000. It'll use more memory, but it might be a little quicker.

[snigdhaAgarwal]

A few questions:

1. Currently the tiff files saved in this manner don't open in several image viewing platforms like ImageJ, giving error like:

libvips saves as a simple flat TIFF by default. Perhaps your viewer is doing nearest neighbour for zoom out?

Try:

dim.write_to_file("x.tif", pyramid=True, tile=True, compression="jpeg")

Deepzoom works well too, if you have a viewer for that.

vipsdisp is a gtk3 viewer that works for most of these formats, though it'll also do nearest neighbour if there's no image pyramid.

https://github.com/jcupitt/vipsdisp

Any ideas on how I can go about debugging this issue?

2. The cropped images are currently in a polygon shape, with the rest of the image being black. Is there a way to make the rest of the area white? I tried fiddling with style and fill in svg viewbox but that didn't work.

[jcupitt]

On 1. I guess this is a JPEG-tiled TIFF, is that right?

Yes, many platforms won't open TIFFs like this, sadly. You'll need to use another format, or another viewer. QuPath works well for this sort of thing, and is a lot more modern than imagej.

[snigdhaAgarwal]

1. Tried opening with Qupath but I get this message and I get a black screen regardless of the choice of Yes or No.
2. Still haven't resolved this. Setting background-color in svg element doesn't work.

[jcupitt]

Have you tried setting pyramid? I see:

$ vips copy summer8.tif x.tif[tile,compression=jpeg,pyramid]
$ vipsheader x.tif
x.tif: 18008x7588 uchar, 3 bands, srgb, tiffload
$ ~/packages/bioformats/QuPath-0.2.3/bin/QuPath-0.2.3 x.tif
...

And it pops up a view window in about a second.

For the background, put the SVG mask into the alpha, then flatten with white.

[jcupitt]

Oh, you're already putting the mask into the alpha, I missed the [3] at the end. Just flatten, eg.:

dim.flatten(background=255).write_to_file('trial2.tiff', pyramid=True, tile=True, compression='jpeg')

That'll flatten out the alpha by multiplying it with white (255).

[jcupitt]

libvips automatically flattens RGBA to RGB when you save TIFF with JPEG compression, so this will work too:

dim.write_to_file('trial2.tiff', pyramid=True, tile=True, compression='jpeg', background=255)