How to get the mse of each pixel in two images of tha same dimension in pyvips

[SikangSHU]

Hello! As my tiff file is big and can't use the way as follows to get the mse of each pixel in two images. Is there any way that using pyvips to get the mse of each pixel in two images of tha same dimension? Thanks!

def mse(imageA, imageB):
    # the 'Mean Squared Error' 
    err = np.sum((imageA.astype("float") - imageB.astype("float")) ** 2)
    err /= float(imageA.shape[0] * imageA.shape[1])
    return err

[SikangSHU]

Hello! As my tiff file is big and can't use the way as follows to get the mse of each pixel in two images. Is there any way that using pyvips to get the mse of each pixel in two images of tha same dimension? Thanks!

def mse(imageA, imageB):
    # the 'Mean Squared Error' 
    err = np.sum((imageA.astype("float") - imageB.astype("float")) ** 2)
    err /= float(imageA.shape[0] * imageA.shape[1])
    return err

I have tried in this way.

img1 = pyvips.Image.tiffload('1.tif')
img2 = pyvips.Image.tiffload('2.tif')

mse_result = ((img1[0] - img2[0]) ** 2 + (img1[1] - img2[1]) ** 2 + (img1[2] - img2[2]) ** 2) / 3

[jcupitt]

Hi again, try:

a = pyvips.Image.new_from_file(xxx)
b = pyvips.Image.new_from_file(yyy)
mse = ((a - b) ** 2).avg()

[jcupitt]

You can stream this computation, so this would be faster:

a = pyvips.Image.new_from_file("a.tif", access="sequential")
b = pyvips.Image.new_from_file("b.tif", access="sequential")
mse = ((a - b) ** 2).avg()
print(f"mse = {mse}")

[SikangSHU]

OK. Thank you very much！

发自我的iPhone

------------------ Original ------------------ From: John Cupitt @.> Date: Fri,Dec 31,2021 5:53 PM To: libvips/pyvips @.> Cc: T.R. @.>, Author @.> Subject: Re: [libvips/pyvips] How to get the mse of each pixel in two images of tha same dimension in pyvips (Issue #296)

You can stream this computation, so this would be faster: a = pyvips.Image.new_from_file("a.tif", access="sequential") b = pyvips.Image.new_from_file("b.tif", access="sequential") mse = ((a - b) ** 2).avg() print(f"mse = {mse}")

— Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you authored the thread.Message ID: @.***>

[SikangSHU]

Hello! Is there any conditional statement in pyvips that I can use to express what I want to express in the code I write as follows? Thank you!

mse_result1 = ((y - y1) ** 2).avg()
mse_result2 = ((y - y2) ** 2).avg()
mse_result3 = ((y - y3) ** 2).avg()

image1 = [ ]
image2 = [ ]
image3 = [ ]

if (mse_result1 > mse_result2 and mse_result3 > mse_result2) :
    stain_space_2_exp = stain_space_2.exp()

    image1 = stain_space_2_exp[0]
    image2 = stain_space_2_exp[1]
    image3 = stain_space_2_exp[2]

elif (mse_result2 > mse_result1 and mse_result3 > mse_result1) :
    stain_space_1_exp = stain_space_1.exp()

    image1 = stain_space_1_exp[0]
    image2 = stain_space_1_exp[1]
    image3 = stain_space_1_exp[2]

else:
    stain_space_3_exp = stain_space_3.exp()

    image1 = stain_space_3_exp[0]
    image2 = stain_space_3_exp[1]
    image3 = stain_space_3_exp[2]

[jcupitt]

Hi again, do you want to do this for each pixel? Or for whole images?

[SikangSHU]

Hi again, do you want to do this for each pixel? Or for whole images?

Hi! I'd like to do this for each pixel of a whole slide image(tiff file format) to get three new images(image1, image2 and image 3).

[jcupitt]

You can use condition.ifthenelse(then_image, else_image) to use a condition image to pick pixels between two images.

You'll need to remove the .avg() from your MSE calculation, since that will find a single number, not a whole image.

For example, to pick the pixel with the lower MSE at each point:

a_mse = ((a - y) ** 2).bandmean()
b_mse = ((b - y) ** 2).bandmean()

out = (a_mse < b_mse).ifthenelse(a, b)

[SikangSHU]

You can use condition.ifthenelse(then_image, else_image) to use a condition image to pick pixels between two images.

You'll need to remove the .avg() from your MSE calculation, since that will find a single number, not a whole image.

For example, to pick the pixel with the lower MSE at each point:

a_mse = ((a - y) ** 2).bandmean()
b_mse = ((b - y) ** 2).bandmean()

out = (a_mse < b_mse).ifthenelse(a, b)

Thanks! I have understood. But how can I pick the pixel if there are five MSE and I'd like to pick the pixel with the lowest MSE at each point. If I write in the way as follows, I can't compare the other four cases to decide which one is the lowest MSE.

mse_result1 = ((y1 - y) ** 2).bandmean()
mse_result2 = ((y2 - y) ** 2).bandmean()
mse_result3 = ((y3 - y) ** 2).bandmean()
mse_result4 = ((y4 - y) ** 2).bandmean()
mse_result5 = ((y5 - y) ** 2).bandmean()

image1  = (mse_result1 > mse_result2 and mse_result3 > mse_result2 and mse_result4 > mse_result2 and mse_result5 > mse_result2).ifthenelse(stain_space_2_exp[0], b)
image2  = (mse_result1 > mse_result2 and mse_result3 > mse_result2 and mse_result4 > mse_result2 and mse_result5 > mse_result2).ifthenelse(stain_space_2_exp[1], b)
image3  = (mse_result1 > mse_result2 and mse_result3 > mse_result2 and mse_result4 > mse_result2 and mse_result5 > mse_result2).ifthenelse(stain_space_2_exp[2], b)

[jcupitt]

You can use the boolean operators to combine conditions, or you can nest ifthenelse.

out = (a < b & a < c).ifthenelse(a, (b < c).ifthenelse(b, c))

Or maybe:

out = (a < b).ifthenelse((a < c).ifthenelse(a, c), (b < c).ifthenelse(b, c))

[SikangSHU]

Thank you! I'll have a try.

[SikangSHU]

Hello! After I use ifthenelse, the type of the image that I get is <bound method Image.ifthenelse of <pyvips.Image 61659x65001 uchar, 1 bands, srgb>> I'd like to do some simple operations on the pixels of the image and convert it to grey16, but an error occurred(TypeError: unsupported operand type(s) for *: 'method' and 'int)when I wrote as follows.

filename = f"image.tif[tile,pyramid,subifd]"
print(f"writing {filename} ...")
image = (255 - image * 255).cast("int").colourspace("grey16")
image.set_type(pyvips.GValue.gstr_type, "image-description",
f"""
......
""")
image.tiffsave('image.tif', tile=True, pyramid=True, subifd=True)

Is there a specific method to output images？Thank you!

[jcupitt]

Something has gone wrong --- you'll get an image back from ifthenelse, not a bound method. Did you forget the ( )?

[SikangSHU]

Something has gone wrong --- you'll get an image back from ifthenelse, not a bound method. Did you forget the ( )?

Hi, I have checked my program and maybe I haven't forgot it.

mse_result1 = ((y1 - y) ** 2).bandmean()
mse_result2 = ((y2 - y) ** 2).bandmean()
mse_result3 = ((y3 - y) ** 2).bandmean()
mse_result4 = ((y4 - y) ** 2).bandmean()
mse_result5 = ((y5 - y) ** 2).bandmean()

gray_CD8 = (mse_result1 < mse_result2 & mse_result1 < mse_result3 & mse_result1 < mse_result4 & mse_result1 < mse_result5).ifthenelse(stain_space_1_exp[0], (mse_result3 < mse_result1 & mse_result3 < mse_result2 & mse_result3 < mse_result4 & mse_result3 < mse_result5).ifthenelse(stain_space_3_exp[0], 1))

filename = f"image.tif[tile,pyramid,subifd]"
print(f"writing {filename} ...")
image = (255 - image * 255).cast("int").colourspace("grey16")
image.set_type(pyvips.GValue.gstr_type, "image-description",
f"""
......
""")
image.tiffsave('image.tif', tile=True, pyramid=True, subifd=True)

The type of mse_result1 is <pyvips.Image 61659x65001 float, 1 bands, srgb> , while that of gray_CD8 is <bound method Image.ifthenelse of <pyvips.Image 61659x65001 uchar, 1 bands, srgb>>.

[jcupitt]

You need to post a complete program I can run.

Here's what you should see:

$ python3
Python 3.9.7 (default, Sep 10 2021, 14:59:43) 
[GCC 11.2.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import pyvips
>>> a = pyvips.Image.new_from_file("k2.jpg")
>>> (a > 128).ifthenelse(128, a)
<pyvips.Image 1450x2048 uchar, 3 bands, srgb>

[SikangSHU]

You need to post a complete program I can run.

Here's what you should see:

$ python3
Python 3.9.7 (default, Sep 10 2021, 14:59:43) 
[GCC 11.2.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import pyvips
>>> a = pyvips.Image.new_from_file("k2.jpg")
>>> (a > 128).ifthenelse(128, a)
<pyvips.Image 1450x2048 uchar, 3 bands, srgb>

OK. the program as follows is a complete program. I know it may not be in line with your idea, I will improve my program after it runs successfully. Thank you!

from numpy import double, linalg
import pyvips
import gc

image = pyvips.Image.tiffload('E:\\testpicture_jupyter\\4colour_ceshitu\\unc9.ome.tif')

image /= 255

y = image.log()

stain_vectors_1 = [
    [0.571, 0.095, 0.767],
    [0.584, 0.258, 0.576],
    [0.577, 0.961, 0.284]
]
stain_inverse_1 = linalg.inv(stain_vectors_1).tolist()
stain_space_1 = y.recomb(stain_inverse_1)
y1 = stain_space_1.recomb(stain_vectors_1)

stain_vectors_2 = [
    [0.095, 0.105, 0.767],
    [0.258, 0.758, 0.576],
    [0.961, 0.644, 0.284]
]
stain_inverse_2 = linalg.inv(stain_vectors_2).tolist()
stain_space_2 = y.recomb(stain_inverse_2)
y2 = stain_space_2.recomb(stain_vectors_2)

stain_vectors_3 = [
    [0.571, 0.767, -0.48],
    [0.584, 0.576, 0.808],
    [0.577, 0.284, -0.343]
]
stain_inverse_3 = linalg.inv(stain_vectors_3).tolist()
stain_space_3 = y.recomb(stain_inverse_3)
y3 = stain_space_3.recomb(stain_vectors_3)

stain_vectors_4 = [
    [0.095, 0.767, -0.553],
    [0.258, 0.576, 0.817],
    [0.961, 0.284, -0.165]
]
stain_inverse_4 = linalg.inv(stain_vectors_4).tolist()
stain_space_4 = y.recomb(stain_inverse_4)
y4 = stain_space_4.recomb(stain_vectors_4)

stain_vectors_5 = [
    [0.105, 0.767, -0.218],
    [0.758, 0.576, 0.649],
    [0.644, 0.284, -0.729]
]
stain_inverse_5 = linalg.inv(stain_vectors_5).tolist()
stain_space_5 = y.recomb(stain_inverse_5)
y5 = stain_space_5.recomb(stain_vectors_5)

stain_space_1_exp = stain_space_1.exp()
stain_space_2_exp = stain_space_2.exp()
stain_space_3_exp = stain_space_3.exp()
stain_space_4_exp = stain_space_4.exp()
stain_space_5_exp = stain_space_5.exp()

mse_result1 = ((y1 - y) ** 2).bandmean()
mse_result2 = ((y2 - y) ** 2).bandmean()
mse_result3 = ((y3 - y) ** 2).bandmean()
mse_result4 = ((y4 - y) ** 2).bandmean()
mse_result5 = ((y5 - y) ** 2).bandmean()

gray_CD8 = []
gray_PanCK = []
gray_PDL1 = []
gray_Hema = []

gray_PanCK = (mse_result2 < mse_result1 & mse_result2 < mse_result3 & mse_result2 < mse_result4 & mse_result2 < mse_result5).ifthenelse
(stain_space_2_exp[0], (mse_result1 < mse_result2 & mse_result1 < mse_result3 & mse_result1 < mse_result4 & mse_result1 < mse_result5).ifthenelse(stain_space_1_exp[1],
(mse_result4 < mse_result1 & mse_result4 < mse_result2 & mse_result4 < mse_result3 & mse_result4 < mse_result5).ifthenelse(stain_space_4_exp[0], 1)))

gray_CD8 = (mse_result1 < mse_result2 & mse_result1 < mse_result3 & mse_result1 < mse_result4 & mse_result1 < mse_result5).ifthenelse
(stain_space_1_exp[0], (mse_result3 < mse_result1 & mse_result3 < mse_result2 & mse_result3 < mse_result4 & mse_result3 < mse_result5).ifthenelse(stain_space_3_exp[0], 1))

gray_PDL1 = (mse_result2 < mse_result1 & mse_result2 < mse_result3 & mse_result2 < mse_result4 & mse_result2 < mse_result5).ifthenelse
(stain_space_2_exp[1], (mse_result5 < mse_result1 & mse_result5 < mse_result2 & mse_result5 < mse_result3 & mse_result5 < mse_result4).ifthenelse(stain_space_5_exp[0], 1))

gray_Hema = (mse_result2 < mse_result1 & mse_result2 < mse_result3 & mse_result2 < mse_result4 & mse_result2 < mse_result5).ifthenelse
(stain_space_2_exp[2], (mse_result1 < mse_result2 & mse_result1 < mse_result3 & mse_result1 < mse_result4 & mse_result1 < mse_result5).ifthenelse(stain_space_1_exp[2],
(mse_result3 < mse_result1 & mse_result3 < mse_result2 & mse_result3 < mse_result4 & mse_result3 < mse_result5).ifthenelse(stain_space_3_exp[1], (mse_result4 < mse_result1 & mse_result4 < mse_result2 &
mse_result4 < mse_result3 & mse_result4 < mse_result5).ifthenelse(stain_space_4_exp[1], stain_space_5_exp[1]))))

filename = f"CD8.tif[tile,pyramid,subifd]"
print(f"writing {filename} ...")
CD8 = (255 - gray_CD8 * 255).cast("int").colourspace("grey16")
CD8.set_type(pyvips.GValue.gstr_type, "image-description",
f"""
<OME xmlns="http://www.openmicroscopy.org/Schemas/OME/2016-06"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
    Creator="OME Bio-Formats 6.7.0"
    xsi:schemaLocation="http://www.openmicroscopy.org/Schemas/OME/2016-06
        http://www.openmicroscopy.org/Schemas/OME/2016-06/ome.xsd">
    <Image ID="Image:0">
        <!-- Minimum required fields about image dimensions -->
        <Pixels DimensionOrder="XYCZT"
                ID="Pixels:0"
                Interleaved="false"
                SizeC="1"
                SizeT="1"
                SizeX="61659"
                SizeY="65001"
                SizeZ="1"
                Type="uint16">
        <Channel Color="-16718848" ID="Channel:0:0" Name="Channel 1" SamplesPerPixel="1"><LightPath/></Channel><TiffData FirstC="0" FirstT="0" FirstZ="0" IFD="0" PlaneCount="1">
        </TiffData>
        </Pixels>
    </Image>
</OME>
""")
CD8.tiffsave('CD8.tif', tile=True, pyramid=True, subifd=True)

filename = f"PanCK.tif[tile,pyramid,subifd]"
print(f"writing {filename} ...")
PanCK = (255 - gray_PanCK * 255).cast("int").colourspace("grey16")
PanCK.set_type(pyvips.GValue.gstr_type, "image-description",
f"""
<OME xmlns="http://www.openmicroscopy.org/Schemas/OME/2016-06"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
    Creator="OME Bio-Formats 6.7.0"
    xsi:schemaLocation="http://www.openmicroscopy.org/Schemas/OME/2016-06
        http://www.openmicroscopy.org/Schemas/OME/2016-06/ome.xsd">
    <Image ID="Image:0">
        <!-- Minimum required fields about image dimensions -->
        <Pixels DimensionOrder="XYCZT"
                ID="Pixels:0"
                Interleaved="false"
                SizeC="1"
                SizeT="1"
                SizeX="61659"
                SizeY="65001"
                SizeZ="1"
                Type="uint16">
        <Channel Color="-16718848" ID="Channel:0:0" Name="Channel 1" SamplesPerPixel="1"><LightPath/></Channel><TiffData FirstC="0" FirstT="0" FirstZ="0" IFD="0" PlaneCount="1">
        </TiffData>
        </Pixels>
    </Image>
</OME>
""")
PanCK.tiffsave('PanCK.tif', tile=True, pyramid=True, subifd=True)

filename = f"PDL1.tif[tile,pyramid,subifd]"
print(f"writing {filename} ...")
PDL1 = (255 - gray_PDL1 * 255).cast("int").colourspace("grey16")
PDL1.set_type(pyvips.GValue.gstr_type, "image-description",
f"""
<OME xmlns="http://www.openmicroscopy.org/Schemas/OME/2016-06"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
    Creator="OME Bio-Formats 6.7.0"
    xsi:schemaLocation="http://www.openmicroscopy.org/Schemas/OME/2016-06
        http://www.openmicroscopy.org/Schemas/OME/2016-06/ome.xsd">
    <Image ID="Image:0">
        <!-- Minimum required fields about image dimensions -->
        <Pixels DimensionOrder="XYCZT"
                ID="Pixels:0"
                Interleaved="false"
                SizeC="1"
                SizeT="1"
                SizeX="61659"
                SizeY="65001"
                SizeZ="1"
                Type="uint16">
        <Channel Color="-16718848" ID="Channel:0:0" Name="Channel 1" SamplesPerPixel="1"><LightPath/></Channel><TiffData FirstC="0" FirstT="0" FirstZ="0" IFD="0" PlaneCount="1">
        </TiffData>
        </Pixels>
    </Image>
</OME>
""")
PDL1.tiffsave('PDL1.tif', tile=True, pyramid=True, subifd=True)

filename = f"Hema.tif[tile,pyramid,subifd]"
print(f"writing {filename} ...")
Hema = (255 - gray_Hema * 255).cast("int").colourspace("grey16")
Hema.set_type(pyvips.GValue.gstr_type, "image-description",
f"""
<OME xmlns="http://www.openmicroscopy.org/Schemas/OME/2016-06"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
    Creator="OME Bio-Formats 6.7.0"
    xsi:schemaLocation="http://www.openmicroscopy.org/Schemas/OME/2016-06
        http://www.openmicroscopy.org/Schemas/OME/2016-06/ome.xsd">
    <Image ID="Image:0">
        <!-- Minimum required fields about image dimensions -->
        <Pixels DimensionOrder="XYCZT"
                ID="Pixels:0"
                Interleaved="false"
                SizeC="1"
                SizeT="1"
                SizeX="61659"
                SizeY="65001"
                SizeZ="1"
                Type="uint16">
        <Channel Color="-16718848" ID="Channel:0:0" Name="Channel 1" SamplesPerPixel="1"><LightPath/></Channel><TiffData FirstC="0" FirstT="0" FirstZ="0" IFD="0" PlaneCount="1">
        </TiffData>
        </Pixels>
    </Image>
</OME>
""")
CD8.tiffsave('Hema.tif', tile=True, pyramid=True, subifd=True)

gc.collect()

[SikangSHU]

I'd like to get four one-band images and add colours on them.