CP3.0: Spline smoothing in CorrectIlluminationCalculate: possibly wrong last row and column


#1

Hello,

i’m building a background illumination correction pipeline and i have found that using splines smoothing works best with my data sets.

(here is an example of the illumination function)

However upon closer inspection, the calculated illumination function shows strange values in the last row and column.

(same example but zoomed on the lower right corner)

These pixels have clearly discrepant values from the rest of the distribution when looking at the histograms.
and of course, substracting the illumination function to the original image brings erroneous pixels in the final image:

So, is it a bug or a property of the spline interpolation (since it appears only with the spline smoothing)?

Thanks.

(here are the pipeline and the original image)

background_corr_32bits.cpproj (615.8 KB)


#2

We found an issue in the “CorrectIluminationCalculate” module from CP2.2 and it appears to be the same in CP3.0 when Rescale the illumination function is set to yes. The CIC module rescale option appears to use a trimmed distribution which can be found at line 907 in the link listed below (pixel_data*robust_factor=0.02). This causes values to be above 1 if rescaling is performed within the CIC module. Our solution was to turn off the rescale option in the CIC module and use a subsequent RescaleIntensity module in the next step which does not use the trimmed distribution, so, the rescaled values do not go above 1 now. Not sure if this is the case with your pipeline since I have not seen it but sound like a similar issue.

The other question that I have wrt your illumination correction is that the shape of your illumination correction image looks to be a little strange. I can see the spinning disk pattern in the image so it appears to not be smoothed much which can increase the pixel s.d. when you divide by the noisy illumination function. The laser illumination should generate a hotspot in the image FOV which looks like a gaussian distribution. We generate our illumination correction images by imaging a conc. dye solution in an adjacent well so the laser illumination pattern is measured for every plate wavelength channel being used on a given day since bumping the fiber could change the laser alignment. I attached a tile scan image of our acquisition from a 638nm dye well to demonstrate what the raw data looks like for generating the 638 illumination correction function. In general, we average ~50 dye images then the illumination correction image is generated from the averaged images on a per channel basis.


#3

Thanks for your answer derek,

unfortunatly, the Rescale option in the CIC module of my pipeline is set to ‘No’. Rather than a trimmed distribution problem, i think the issue is coming from the call of «affine_transform» from the scipy.ndimage module, in the centrosome.bg_compensate module:

(https://github.com/CellProfiler/centrosome/blob/master/centrosome/bg_compensate.py#L348).

Indeed, from the doc of affine_transform, the keywords «mode» and «cval» will pad the results with zeros when the output coordinates points outside the boundaries of the input:

extract from scipy,ndimage.affine_transform docs
«mode : str, optional
Points outside the boundaries of the input are filled according to the given mode (‘constant’, ‘nearest’, ‘reflect’, ‘mirror’ or ‘wrap’). Default is ‘constant’.
cval : scalar, optional
Value used for points outside the boundaries of the input if mode=‘constant’. Default is 0.0»

And upon testing the backgr function (from centrosome.bg_compensate) with the test image i joined in my original post, and with the settings form my pipeline (all the keywords of backgr are the defaults ones except for “scale=2”), i indeed found the last column and the last row of the output set to 0.

here is the last part of the backgr function which i slightly modified in order to check the indexing and shapes:

…:
output = np.zeros(orig_shape, img.dtype)
print(‘RES SHAPE’,res.shape)
print(‘output_shape =’, tuple(clip_shape))
print(‘clips’,clip_imin,clip_imax,’,’, clip_jmin,clip_jmax)

aff = affine_transform(res, inverse_transform, 
                       output_shape = tuple(clip_shape),
                       order = 3)
print('affshape',aff.shape)
plt.imshow(aff[-10:,-10:])
output[clip_imin:clip_imax, clip_jmin:clip_jmax] = aff
if input_mask is not None:
    output[~input_mask] = 0
return output

and here is the call of the function and the results:

In [135]: outp = backgr(img, scale=2)
RES SHAPE (375, 375)
output_shape = (750, 750)
clips 0 750 , 0 750
affshape (750, 750)

image

(the yellow pixels are all equal to 0, the others are negative hee)

But i don’t really get why this is happenning, since as you can see the indexing is correct, and when manually calculating the correspondance between output and input coordinates, as explained in the affine_transform doc:

«Given an output image pixel index vector o, the pixel value is determined from the input image at position np.dot(matrix,o) + offset.»

everything seems fine (the offset in the backgr function is null):

In [136]: a=np.array([749,749])
In [137]: np.dot(inverse_transform,a)
array([ 374., 374.])

where inverse_transform is

In [138]: inverse_transform
array([[ 0.49933244, 0. ],
[ 0. , 0.49933244]])

Any idea what is going on? Could it be a floating-point rounding issue or something like it?

As for your others questions, i am not working with intensity images, but with phase images reconstructed from a digital hologram without any dye. This explains why i don’t have a regular illumination function since it is not really an illumination function, and why there are concentric patterns which are interference patterns. But the principle of the background correction still remains valid/usefull.

Finally i also have a bunch of questions/remarks related to cellprofiler. I’ll ask them here but will start a new thread if needed:

  1. Why can’t we save, with the saveImage module, 32 bits images with pixels values outside the range of [-1,1]? (i get the neccessity to rescale for the treatments but why for the save as well?)

  2. When i compare the given WALL_TIME per worker and per module in the console, with the actual time it takes to treat one set of data, i get a real time twice as long as the given wall time. What i am missing here?

  3. I tried to use cellProfiler as a python package and was wondering if there was an available documentation to do so (the wiki page: https://github.com/CellProfiler/CellProfiler/wiki/CellProfiler-as-a-Python-package is a bit short )?

  4. same question for headless mode.

thanks for your time and great work.


#4

Hi,

That does indeed look weird, I expect it’s a bug. Would you mind making a GitHub issue?


#5

Hi bcimini,
It’s done


#7

Oh that’s why your illumination correction images look different than what I expected.