In [1]:
import diplib as dip
import random
import numpy as np
This function creates two test images composed of blobs, with a fraction overlap of the dots overlapping.
In [2]:
def generate_images(overlap):
sd = 0.001 # noise std. dev.
sz = 5.0 # size of dot (sigma of Gaussian)
scale = sz*sz*2*3.14159
channel1 = dip.Image([256,256], 1, 'SFLOAT')
channel2 = dip.Image([256,256], 1, 'SFLOAT')
channel1.Fill(0)
channel2.Fill(0)
for jj in range(8):
for ii in range(8):
x = 4*sz + ii * 6*sz
y = 4*sz + jj * 6*sz
dip.DrawBandlimitedPoint(channel1, [x, y], scale, sz)
if ii < 6:
# If larger, the 2nd channel doesn't have a dot
if jj * 8 + ii > 8 * 8 * overlap:
# Not overlapping points, move them!
x += 3*sz
y += 1*sz
dip.DrawBandlimitedPoint(channel2, [x, y], [scale], [sz])
channel1 = dip.ClipLow(dip.GaussianNoise(channel1, sd**2), 0)
channel2 = dip.ClipLow(dip.GaussianNoise(channel2, sd**2), 0)
return channel1, channel2
Let's test see what these images look like:
In [3]:
channel1, channel2 = generate_images(0.5)
dip.JoinChannels([channel1, channel2]).Show()
Now we'll run through various values of overlap, generate images, and compute colocalization:
In [4]:
for overlap in [0.1, 0.3, 0.5, 0.7, 0.9]:
channel1, channel2 = generate_images(overlap)
print()
print(overlap)
print('PearsonCorrelation: ', round(dip.PearsonCorrelation(channel1, channel2), 3))
print('MandersOverlapCoefficient: ', round(dip.MandersOverlapCoefficient(channel1, channel2), 3))
print('IntensityCorrelationQuotient: ', round(dip.IntensityCorrelationQuotient(channel1, channel2), 3))
coef = dip.MandersColocalizationCoefficients(channel1, channel2, None, 0.2, 0.2)
print('MandersColocalizationCoefficients: ', round(coef[0], 3), round(coef[1], 3))
coef = dip.CostesColocalizationCoefficients(channel1, channel2)
print('CostesColocalizationCoefficients: ', round(coef[0], 3), round(coef[1], 3))
To apply Costes' significance test, we need to estimate the width of the autocorrelation function
In [5]:
ac = dip.AutoCorrelationFT(channel1)
ac = ac > dip.Maximum(ac)[0][0] / 2 # half maximum
ac = dip.Label(ac)
cc = dip.GetImageChainCodes(ac, ac[ac.Size(0)//2, ac.Size(1)//2]) # find central blob
bb = cc[0].BoundingBox()
blockSizes = [bb[1][0] - bb[0][0], bb[1][1] - bb[0][1]] # width and height, corresponds to full width at half maximum
print(blockSizes)
This width we use as the blockSizes parameter to dip.CostesSignificaneTest. We explore the region of the overlap values that are interesting:
In [6]:
for overlap in [0.00, 0.14, 0.16, 0.18, 0.20, 0.22]:
channel1, channel2 = generate_images(overlap)
print(overlap, round(dip.CostesSignificanceTest(channel1, channel2, None, blockSizes, 500), 4))
In [ ]: