libvips / pyvips

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

How to stitch the images of whole slide. #368

Open kuldeep203 opened 1 year ago

kuldeep203 commented 1 year ago

Currently i am working on digital microscope. I have to scan whole slide and then stitch them together to create a single image. i am scanning the slide in snake manner as shown in the attachment. Currently i have 19 row and 17 column. Please suggest me how i can achieve that. Screenshot from 2023-01-02 15-13-56

jcupitt commented 1 year ago

Hi @kuldeep203,

I would use python and either assemble with join or mosaic, depending on the accuracy of your microscope stage.

You will probably also need to do radiometric (lighting uniformity, lens vignetting) correction, and possibly geometric (lens distortion).

kuldeep203 commented 1 year ago

can you tell me how i join arrayjoin and mosaic for row and column wise. take the above example of screen shot. how i can define overlap in array join method.

jcupitt commented 1 year ago

Use the hspacing and vspacing arguments to arrayjoin to set the gap between tiles.

For example, make a set of test data with:

vips dzsave nina.jpg x --depth one --tile-width 400 --tile-height 400 --overlap 50

Now centre tiles are 500x500, with 100 pixels of overlap along every edge. Tiles are named as {x}_{y}.jpeg. Reassemble with:

cd x_files/0
vips arrayjoin "$(ls *.jpeg | sort -t_ -k2g -k1g)" x.jpg --across 16 --hspacing 400 --vspacing 400

This is pretty ugly with bash, you'd find python is simpler.

kuldeep203 commented 1 year ago

hi @jcupitt Thanks for the reply, i used arrayjoin but i am getting grid like structure as shown in the attachment. How i can color balance the image after stitching. Screenshot from 2023-01-04 09-43-29

jcupitt commented 1 year ago

That's lens vignetting, so you need to flatfield each capture.

Take a pic of just the white background (no slide), then calculate something like:

blur = white.gaussblur(10)
flatfield = (blur.max() / blur).copy_memory()

Now flatfield will be 1.0 for the brightest part of the field, and maybe 1.4 at the edges. Flatten each capture with:

frame = (frame * flatfield).cast("uchar")

Now join the frames and it should look better.

jcupitt commented 1 year ago

To get a good result you will need to consider and correct other sources of error, such as camera linearity, cast correction, black levels, changes in illumination strength, optical scattering, dead pixels, errors in stage movement, optical distortion, errors in stage and camera alignment, etc.

kuldeep203 commented 1 year ago

Thanks for the reply. i have some question. Is gaussblur is in inbuilt in pyvips or you are using this from opencv. what does frame refer too. is it refer to images. image = frame. images = frames. and white = white image captured.

jcupitt commented 1 year ago

Yes, that's all pyvips, and everything is an image.

kuldeep203 commented 1 year ago

vips dzsave nina.jpg x --depth one --tile-width 400 --tile-height 400 --overlap 50

With the above command center tile will be 450 or 500 because i am getting 450X450 not as you mentioned 500X500.

and second thing is that i am not able to calculate exact overlap between two captured images. Is there anyway to overcome this like we can give approximate overlap if cant then how we can calculate exact overlap between two images. Sorry to ask too many question as i am new to image processing. Thanks for your continue help.

jcupitt commented 1 year ago

Yes, the edge tiles will be 450x450, the centre tiles will be 500x500.

arrayjoin only does exact grids. If you need to correct for motion errors, you can use the mosaic functions to build the grid as a set of left-right and top-bottom joins. These will search for the exact overlap for you, then do a feathered join.

They will not work well for very large grids (more than 20x20?), since errors will tend to accumulate. In this case, you will need to assemble in two passes: first search all overlaps for the exact offset, then make a big linear system with the constraints, then solve using least mean square to find the set of tile positions which minimise the average offset error, and finally do a pair-wise assembly using the computed positions.

jcupitt commented 1 year ago

... a lot depends on how much control you have over the microscope. Can you adjust the camera pitch, yaw and roll, for example? Can you adjust the stage as well? How well do images align if you swap the objective?

kuldeep203 commented 1 year ago

In this case, you will need to assemble in two passes: first search all overlaps for the exact offset, then make a big linear system with the constraints, then solve using least mean square to find the set of tile positions which minimise the average offset error, and finally do a pair-wise assembly using the computed positions. Can you provide me the example of it.Because i am new at image processing. i don't know too much about it. And second thing is that at last i will have approx 5000 images. So i can get exact mosaic if you want i can share you some input images.

jcupitt commented 1 year ago

Sorry, I don't have time to do the work for you. Find a local expert to help with the maths.

kuldeep203 commented 1 year ago

Thanks @jcupitt

kuldeep203 commented 1 year ago

hi @jcupitt Please look at the pic and let me know i can remove visible lines between the pic. i applied the flat field on every images as well as the gauss blur filter. i use merge operation for making the mosaic. x38 Please find the code and the attachment of image. `#!/usr/bin/env python

import sys import pyvips from datetime import datetime

now = datetime.now()

current_time = now.strftime("%H:%M:%S") print("Current Time =", current_time)

overlap joins by this many pixels

H_OVERLAP = 200 V_OVERLAP = 100

number of images in mosaic

ACROSS = 13 DOWN = 17

if len(sys.argv) < 2 + ACROSS * DOWN: print('usage: %s output-image input1 input2 ..') sys.exit(1)

def join_left_right(filenames): images = [pyvips.Image.new_from_file(filename) for filename in filenames] row = images[0] for image in images[1:]: image.addalpha() row = row.merge(image, 'horizontal', H_OVERLAP - row.width, 0,mblend = 100)

row.globalbalance()

return row

def join_top_bottom(rows): image = rows[0] for row in rows[1:]: print(row) image = image.merge(row, 'vertical', 0, V_OVERLAP - image.height,mblend = 400)

image.globalbalance()

return image

rows = [] for y in range(0, DOWN):

print(y)

start = 2 + y * ACROSS

end = start + ACROSS
images_list = sys.argv[start:end]

if y%2!=0:
    print(y)
    print(images_list)
    images_list1 = images_list.reverse()
    print("else")
    print(images_list1)
    rows.append(join_left_right(images_list))
else:
    rows.append(join_left_right(images_list))

image = join_top_bottom(rows) image.globalbalance() image.write_to_file(sys.argv[1])

`

jcupitt commented 1 year ago

Hello, you've not flatfielded the frames. Post some sample data and I'll try for you.

kuldeep203 commented 1 year ago

I have flat fielded the frames that is different file. i flat field all images and write to another folder and then use that folder for stitching. ` def list_files(file_filder): arr = []

folder path

dir_path = file_filder
# list to store files
# Iterate directory
for path in os.listdir(file_filder):
    # check if current path is a file
    if os.path.isfile(os.path.join(dir_path, path)):
        arr.append(dir_path + path)
return arr

def make_image_flat_filed(): tiff_files_li = list_files("/home/kuldeep/Downloads/Inferencing_with_Button/Input/") tiff_files_li.sort() f = image_correct() for i in tiff_files_li: k = i.replace("/home/kuldeep/Downloads/Inferencing_with_Button/Input/","") print(i) tile = pyvips.Image.new_from_file(i) img = (tile * f).cast("uchar") img.write_to_file("/home/kuldeep/motherson/priid-python/Output/{}".format(k))

` Please find link of google drive for images. Image data are 17X13 (17 row and 13 column).

jcupitt commented 1 year ago

Your tiles have not been flatfielded, something has gone wrong in that part of the process. Share your white image plus one tile of the slide, direct from the camera.

kuldeep203 commented 1 year ago

Please find the white image IMG202301051157362 as well as the input image. image_009 The white image is taken by mobile phone as if am taking the white paper pic i am not getting white image.Please find the attach image. MicrosoftTeams-image

jcupitt commented 1 year ago

You need to take the white image using the microscope, using the same lens, focus and aperture as you use for scanning.

I would do it as part of slide acquisition: put the slide on the microscope, focus it and set the exposure, then drive the stage to the side somewhere so you have empty glass, and take a pic. That image can then be used to correct for lighting and lens effects for your sample frames.

jcupitt commented 1 year ago

Your input image looks out of focus. It's sharp in the bottom middle, but very blurry at the edges.

You will need to fix this as well. Perhaps you can find a better lens, or maybe reduce the field of view to just the area of good pixels.

kuldeep203 commented 1 year ago

Hi @jcupitt i have corrected the lenses and the quality of the images is good and merging function is also good. Now i am trying to stitch 10000 image with merge function but not able to do. My program is stuck at one place and i not getting any memory issue. So please let me know how i can improve that.

jcupitt commented 1 year ago

You'll need to explain how you are merging them together. Post some demo code that illustrates the problem.

kuldeep203 commented 1 year ago

Please find the code below:-


import sys
import pyvips
from datetime import datetime

now = datetime.now()

current_time = now.strftime("%H:%M:%S")
print("Current Time =", current_time)
# overlap joins by this many pixels
H_OVERLAP = 100

V_OVERLAP = 100

# number of images in mosaic
ACROSS = 100
DOWN = 80

if len(sys.argv) < 2 + ACROSS * DOWN:
    print('usage: %s output-image input1 input2 ..')
    sys.exit(1)

def join_left_right(filenames):
    images = [pyvips.Image.new_from_file(filename) for filename in filenames]
    row = images[0]
    for image in images[1:]:
        image.addalpha()
        row = row.merge(image, 'horizontal', H_OVERLAP - row.width, 0, mblend=20)
        # row.globalbalance()
    return row

def join_top_bottom(rows):
    image = rows[0]
    for row in rows[1:]:
        print(row)
        image = image.merge(row, 'vertical', 0, V_OVERLAP - image.height, mblend=20)
        # image.globalbalance()
    return image

rows = []
for y in range(0, DOWN):

    start = 2 + y * ACROSS
    end = start + ACROSS
    # print("start is {} and and is {} is =".format(start, end))
    images_list = sys.argv[start:end]
    # images_list.reverse()

    if y % 2 != 0:

        images_list1 = images_list.reverse()

        rows.append(join_left_right(images_list))
    else:
        rows.append(join_left_right(images_list))

image = join_top_bottom(rows)
image.globalbalance()
# image.gaussblur(5)
image.write_to_file(sys.argv[1], tile=True, compression="jpeg")

**`**

system info: Architecture: x86_64 CPU op-mode(s): 32-bit, 64-bit Address sizes: 39 bits physical, 48 bits virtual Byte Order: Little Endian CPU(s): 8 On-line CPU(s) list: 0-7 Vendor ID: GenuineIntel Model name: Intel(R) Core(TM) i5-10300H CPU @ 2.50GHz CPU family: 6 Model: 165 Thread(s) per core: 2 Core(s) per socket: 4 Socket(s): 1 Stepping: 2 CPU max MHz: 4500.0000 CPU min MHz: 800.0000 BogoMIPS: 4999.90

Image Size: = 896,684

Thanks in advance for help.

kuldeep203 commented 1 year ago

In this case, you will need to assemble in two passes: first search all overlaps for the exact offset, then make a big linear system with the constraints, then solve using least mean square to find the set of tile positions which minimise the average offset error, and finally do a pair-wise assembly using the computed positions.

Hi i am trying to find the exact overlap with the code mention below

#!/usr/bin/python

import sys
import pyvips

# you can set the size of the window as well, but 10x10 pixels (the default) is
# probably fine

left = pyvips.Image.new_from_file(sys.argv[1])
right = pyvips.Image.new_from_file(sys.argv[2])

# the expected overlap, or the number of pixels the two images have in common
overlap = int(sys.argv[3])

# harea is the maximum displacement we detect, in pixels ... this is half
# the size of the area searched

# bandno is the band we search in ... green (band 1) is best for JPG images,
# since it has double the resolution

# dx0/dy0 are optional output args ... we set them True and pyvips will return
# their rvalues in result

join, optional = left.mosaic(right, 'horizontal',
                             left.width - overlap, 0, 0, 0,
                             harea=15, bandno=0,
                             dx0=True, dy0=True)

print 'horizontal offset =', optional['dx0']
print 'vertical offset =', optional['dy0']

print 'writing join to x.png ...'
join.write_to_file('x.png')

this is woking fine , but you are talking about least mean square to find the set of tile positions which minimise the average offset error can you provide me the example of that. currently i am stitching 90000 images and able to stitch but due to movement of the stage under camera i am not able to give the exact overlap. so i can resolve this.

jcupitt commented 1 year ago

Sorry, I don't have time to code it up, maybe ask a maths friend? It's an over-determined linear system, so you just need to express the set of discovered overlaps as a big matrix and then invert it.