BioImage slicing and stitching in Python without losing information

Question

I am a python newbie learning by doing and working as a medical researcher on large images with a lot of biological data about cell marker expressions. I have a short processing code that does a channel separation and turning the channels as fluorescene (pseudo) and/or grayscale. When trying to process a larger image we dont have enough memory so the way to go is to slice the image in processable tiles, perform the channel formatting and then stich the image back together. I really have 0 idea how the stitching and slicing work incorporated with this channel separation.

I used these snippets from scikit_image

import numpy as np`
import matplotlib.pyplot as plt

from skimage import data
from skimage.color import rgb2hed, hed2rgb

**# Example IHC image(512, 512, 3) uint8**
ihc_rgb = data.immunohistochemistry()

**# Separate the stains from the IHC image**
ihc_hed = rgb2hed(ihc_rgb)

**# Create an RGB image for each of the stains**
null = np.zeros_like(ihc_hed[:, :, 0])
ihc_h = hed2rgb(np.stack((ihc_hed[:, :, 0], null, null), axis=-1))
ihc_e = hed2rgb(np.stack((null, ihc_hed[:, :, 1], null), axis=-1))
ihc_d = hed2rgb(np.stack((null, null, ihc_hed[:, :, 2]), axis=-1))

**# Display**
fig, axes = plt.subplots(2, 2, figsize=(7, 6), sharex=True, sharey=True)
ax = axes.ravel()

ax[0].imshow(ihc_rgb)
ax[0].set_title("Original image")

ax[1].imshow(ihc_h)
ax[1].set_title("Hematoxylin")

ax[2].imshow(ihc_e)
ax[2].set_title("Eosin")  # Note that there is no Eosin stain in this image

ax[3].imshow(ihc_d)
ax[3].set_title("DAB")

for a in ax.ravel():
    a.axis('off')

fig.tight_layout()

from skimage.exposure import rescale_intensity

**# Rescale hematoxylin and DAB channels and give them a fluorescence look**
h = rescale_intensity(ihc_hed[:, :, 0], out_range=(0, 1),
                      in_range=(0, np.percentile(ihc_hed[:, :, 0], 99)))
d = rescale_intensity(ihc_hed[:, :, 2], out_range=(0, 1),
                      in_range=(0, np.percentile(ihc_hed[:, :, 2], 99)))
**
**# Cast the two channels into an RGB image, as the blue and green channels
**# respectively******
zdh = np.dstack((null, d, h))

fig = plt.figure()
axis = plt.subplot(1, 1, 1, sharex=ax[0], sharey=ax[0])
axis.imshow(zdh)
axis.set_title('Stain-separated image (rescaled)')
axis.axis('off')
plt.show()

I want to translate this to a large image.

Answer 1

given the advice in how to read a large image by lines with python :

I would suggest you tile your rather large TIFF outside of Python before you start. Drawing heavily on John's (@jcupitt) answer here, you can do that with vips.

I tried to reproduce your problem. I am a newbie to image-processing and Python too, so bear with me , my code below, I choose the pyvips that wraps the libvips image processing library because I found a lot of help about it here on SO:

from PIL import Image
import numpy as np

import pyvips

# CREATE 10000*10000*3 black image 

# array = np.zeros((10000, 10000, 3), dtype=np.uint8)
# image = pyvips.Image.new_from_array(array)
# image.write_to_file("huge.tif")


# fields = image.get_fields()
 
# print(fields) 


# img = image.numpy()

# print('img : ', img.size, img.shape , img.ndim)


# load the 10000*10000*3 black image , without loaing it !!! this should be where we can load your huge tif
img_loaded = pyvips.Image.new_from_file('huge.tif', access='sequential')


print(img_loaded, img_loaded.width , img_loaded.height, img_loaded.bands)


hg = 0

for i in range(img_loaded.height//1000):
    
    # https://stackoverflow.com/questions/48281086/extracting-a-region-of-interest-from-an-image-file-without-reading-the-entire-im
    
    
    # uses pyvips pointer to loaded image to load part of it 
    # get the numpy array corresponding to it
    # change same of the images
    # saves them one by one 
    #IMPORTANT script works because we know 10000/10 is round no pixel leftover
    
    a = img_loaded.crop(0, i*1000, img_loaded.width, 1000)
    
    
    arr = a.numpy()
    
    print(arr.size , arr.shape , arr.ndim)
    
    # aa = Image.fromarray(arr)
    # aa.show()
    
    if i%2 == 0 :
        
        arr[:,:,0] = 0
        arr[:,:,1] = 255
        arr[:,:,2] = 0
        
        print('i%2 : ', i % 2)
        
        aa = Image.fromarray(arr)
        aa.show()
        
        aa.save('bit_'+str(i)+'.tif')
        
    
    elif i == 3 or i == 9 :
        
        arr[:,:,0] = 255
        arr[:,:,1] = 0
        arr[:,:,2] = 0
        
    
        aa = Image.fromarray(arr)
        aa.show()
        
        aa.save('bit_'+str(i)+'.tif')
       
    else:
        
        print('i%2 : ', i % 2)
        
        aa = Image.fromarray(arr)
        aa.show()
        
        aa.save('bit_'+str(i)+'.tif')
        
 
        
# https://stackoverflow.com/questions/50297176/merge-large-images-on-disk

# load each of the images from above without loading them , and join again

images = [pyvips.Image.new_from_file(filename, access="sequential")
          for filename in ['bit_'+str(i)+'.tif' for i in range(img_loaded.height//1000)   ] ]

final = pyvips.Image.arrayjoin(images, across = 1)

final.write_to_file('result.tiff')

the first lines of my script:

# CREATE 10000*10000*3 black image 

# array = np.zeros((10000, 10000, 3), dtype=np.uint8)
# image = pyvips.Image.new_from_array(array)
# image.write_to_file("huge.tif")

create my supposed huge.tif file , is not that big but ......

exiftool huge.tif gives the details :

Image Width                     : 10000
Image Height                    : 10000
Bits Per Sample                 : 8 8 8

the second part of the code load , without loading in memory , such file

and operates same changes on part of it extracted and converted to numpy arrays that are finally saved as smaller .tif files via Pillowlibrary.

Then the resulting pics are not loaded again via pyvips and merged to a not in memory file written to disk, I don't know if I made same mistake the resulting result.tiff files has same dimensions of input one :

I am new to big tiff files , this should provide a way to handle them, see inside my code some link to relevant SO questions whose solutions I used to create this snippet