From 936bda0aec75a0e978638da5a9f2b7cea0d0a80b Mon Sep 17 00:00:00 2001 From: George Emanuel Date: Fri, 17 Apr 2020 16:12:30 -0400 Subject: [PATCH 01/12] Incremented version. --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index a45ed54b..25090285 100644 --- a/setup.py +++ b/setup.py @@ -17,7 +17,7 @@ setuptools.setup( name="merlin", - version="0.1.6", + version="0.1.7", description="MERFISH decoding software", author="George Emanuel", author_email="emanuega0@gmail.com", From bb97a6fb616173cd1e1f0ad05ae8f3a7c6657cc5 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Tue, 28 Apr 2020 17:51:28 -0400 Subject: [PATCH 02/12] Add utility functions for calculationg high/lowpass images and for calculating pixel significance from foreground and background image estimates. --- merlin/util/imagefilters.py | 67 +++++++++++++++++++++++++++++++++++-- 1 file changed, 64 insertions(+), 3 deletions(-) diff --git a/merlin/util/imagefilters.py b/merlin/util/imagefilters.py index 15e62185..aa85a4dd 100755 --- a/merlin/util/imagefilters.py +++ b/merlin/util/imagefilters.py @@ -6,6 +6,67 @@ """ +def est_significance(foreground: np.ndarray, + background: np.ndarray) -> np.ndarray: + """ + Args: + foreground: the image foreground estimate. + background: the image background estimate. + + Returns: + an estimate of the significance of each pixel in units of + sigma (or standard deviation). + """ + snbfIm = foreground/np.sqrt(background) + + # Rescale snbfIM so that it's histogram is Gaussian with + # zero offset and sigma = 1.0. + # + + # Why subtract 1.0 in order to make offset = 0.0? This is + # emperical, not entirely sure I understand why the offset + # is 1.0. + # + snbfIm -= 1.0 + + sx = min(snbfIm.shape[0], 512) + sy = min(snbfIm.shape[1], 512) + + [hist, edges] = np.histogram(snbfIm[0:sx, 0:sy], bins=100, range=(0, 6)) + centers = 0.5*(edges[1:] + edges[:-1]) + + sigma = np.sum(hist*centers)/np.sum(hist) + + return snbfIm/sigma + + +def high_low_filter(image: np.ndarray, + windowSize: int, + sigma: float, + reps: int) -> np.ndarray: + """ + Args: + image: the input image to be filtered + windowSize: the size of the Gaussian kernel to use. + sigma: the sigma of the Gaussian. + reps: number of repetitions of foreground/background estimation. + + Returns: + the high pass and low pass filtered images. + """ + lowpass = np.copy(image) + for i in range(reps): + lowpass = cv2.GaussianBlur(lowpass, + (windowSize, windowSize), + sigma, + borderType=cv2.BORDER_REPLICATE) + highpass = image - lowpass + highpass[lowpass > image] = 0 + lowpass = image - highpass + + return [highpass, lowpass] + + def high_pass_filter(image: np.ndarray, windowSize: int, sigma: float) -> np.ndarray: @@ -23,6 +84,6 @@ def high_pass_filter(image: np.ndarray, (windowSize, windowSize), sigma, borderType=cv2.BORDER_REPLICATE) - gauss_highpass = image - lowpass - gauss_highpass[lowpass > image] = 0 - return gauss_highpass + highpass = image - lowpass + highpass[lowpass > image] = 0 + return highpass From ea70af10bdc1a27b5a827ec3fdf25cfaca298140 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Thu, 30 Apr 2020 17:36:45 -0400 Subject: [PATCH 03/12] Add clipping for small / negative background values. Remove sigma scaling as this is not necessary if the gain is set correctly. --- merlin/util/imagefilters.py | 18 +++++------------- 1 file changed, 5 insertions(+), 13 deletions(-) diff --git a/merlin/util/imagefilters.py b/merlin/util/imagefilters.py index aa85a4dd..9e05e789 100755 --- a/merlin/util/imagefilters.py +++ b/merlin/util/imagefilters.py @@ -17,11 +17,11 @@ def est_significance(foreground: np.ndarray, an estimate of the significance of each pixel in units of sigma (or standard deviation). """ - snbfIm = foreground/np.sqrt(background) - - # Rescale snbfIM so that it's histogram is Gaussian with - # zero offset and sigma = 1.0. + # Clip so that we don't try and take the sqrt of 0.0 or + # a negative number. # + background = np.clip(background, 1.0, None) + snbfIm = foreground/np.sqrt(background) # Why subtract 1.0 in order to make offset = 0.0? This is # emperical, not entirely sure I understand why the offset @@ -29,15 +29,7 @@ def est_significance(foreground: np.ndarray, # snbfIm -= 1.0 - sx = min(snbfIm.shape[0], 512) - sy = min(snbfIm.shape[1], 512) - - [hist, edges] = np.histogram(snbfIm[0:sx, 0:sy], bins=100, range=(0, 6)) - centers = 0.5*(edges[1:] + edges[:-1]) - - sigma = np.sum(hist*centers)/np.sum(hist) - - return snbfIm/sigma + return snbfIm def high_low_filter(image: np.ndarray, From 325e96ade440f7c6f2d9d501288453e47319cb71 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Thu, 30 Apr 2020 17:38:06 -0400 Subject: [PATCH 04/12] Add EstimatePixelSignificance class. Move some methods from DeconvolutionPreprocess into Preprocess so that they can also be used by EstimatePixelSignificance. --- merlin/analysis/preprocess.py | 136 +++++++++++++++++++++++++--------- 1 file changed, 99 insertions(+), 37 deletions(-) diff --git a/merlin/analysis/preprocess.py b/merlin/analysis/preprocess.py index 32a904f9..8c718842 100755 --- a/merlin/analysis/preprocess.py +++ b/merlin/analysis/preprocess.py @@ -14,53 +14,22 @@ class Preprocess(analysistask.ParallelAnalysisTask): """ An abstract class for preparing data for barcode calling. """ + + def __init__(self, dataSet, parameters=None, analysisName=None): + super().__init__(dataSet, parameters, analysisName) + + if 'codebook_index' not in self.parameters: + self.parameters['codebook_index'] = 0 def _image_name(self, fov): destPath = self.dataSet.get_analysis_subdirectory( self.analysisName, subdirectory='preprocessed_images') return os.sep.join([destPath, 'fov_' + str(fov) + '.tif']) - def get_pixel_histogram(self, fov=None): - if fov is not None: - return self.dataSet.load_numpy_analysis_result( - 'pixel_histogram', self.analysisName, fov, 'histograms') - - pixelHistogram = np.zeros(self.get_pixel_histogram( - self.dataSet.get_fovs()[0]).shape) - for f in self.dataSet.get_fovs(): - pixelHistogram += self.get_pixel_histogram(f) - - return pixelHistogram - def _save_pixel_histogram(self, histogram, fov): self.dataSet.save_numpy_analysis_result( histogram, 'pixel_histogram', self.analysisName, fov, 'histograms') - -class DeconvolutionPreprocess(Preprocess): - - def __init__(self, dataSet, parameters=None, analysisName=None): - super().__init__(dataSet, parameters, analysisName) - - if 'highpass_sigma' not in self.parameters: - self.parameters['highpass_sigma'] = 3 - if 'decon_sigma' not in self.parameters: - self.parameters['decon_sigma'] = 2 - if 'decon_filter_size' not in self.parameters: - self.parameters['decon_filter_size'] = \ - int(2 * np.ceil(2 * self.parameters['decon_sigma']) + 1) - if 'decon_iterations' not in self.parameters: - self.parameters['decon_iterations'] = 20 - if 'codebook_index' not in self.parameters: - self.parameters['codebook_index'] = 0 - - self._highPassSigma = self.parameters['highpass_sigma'] - self._deconSigma = self.parameters['decon_sigma'] - self._deconIterations = self.parameters['decon_iterations'] - - self.warpTask = self.dataSet.load_analysis_task( - self.parameters['warp_task']) - def fragment_count(self): return len(self.dataSet.get_fovs()) @@ -76,6 +45,18 @@ def get_dependencies(self): def get_codebook(self) -> codebook.Codebook: return self.dataSet.get_codebook(self.parameters['codebook_index']) + def get_pixel_histogram(self, fov=None): + if fov is not None: + return self.dataSet.load_numpy_analysis_result( + 'pixel_histogram', self.analysisName, fov, 'histograms') + + pixelHistogram = np.zeros(self.get_pixel_histogram( + self.dataSet.get_fovs()[0]).shape) + for f in self.dataSet.get_fovs(): + pixelHistogram += self.get_pixel_histogram(f) + + return pixelHistogram + def get_processed_image_set( self, fov, zIndex: int = None, chromaticCorrector: aberration.ChromaticCorrector = None @@ -100,6 +81,32 @@ def get_processed_image( chromaticCorrector) return self._preprocess_image(inputImage) + +class DeconvolutionPreprocess(Preprocess): + + def __init__(self, dataSet, parameters=None, analysisName=None): + super().__init__(dataSet, parameters, analysisName) + + if 'highpass_sigma' not in self.parameters: + self.parameters['highpass_sigma'] = 3 + if 'decon_sigma' not in self.parameters: + self.parameters['decon_sigma'] = 2 + if 'decon_filter_size' not in self.parameters: + self.parameters['decon_filter_size'] = \ + int(2 * np.ceil(2 * self.parameters['decon_sigma']) + 1) + if 'decon_iterations' not in self.parameters: + self.parameters['decon_iterations'] = 20 + if 'codebook_index' not in self.parameters: + self.parameters['codebook_index'] = 0 + + self._highPassSigma = self.parameters['highpass_sigma'] + self._deconSigma = self.parameters['decon_sigma'] + self._deconIterations = self.parameters['decon_iterations'] + + self.warpTask = self.dataSet.load_analysis_task( + self.parameters['warp_task']) + + def _high_pass_filter(self, inputImage: np.ndarray) -> np.ndarray: highPassFilterSize = int(2 * np.ceil(2 * self._highPassSigma) + 1) hpImage = imagefilters.high_pass_filter(inputImage, @@ -162,3 +169,58 @@ def _preprocess_image(self, inputImage: np.ndarray) -> np.ndarray: filteredImage, deconFilterSize, self._deconSigma, self._deconIterations).astype(np.uint16) return deconvolvedImage + + +class EstimatePixelSignificance(Preprocess): + """ + Estimates pixel significance in units of sigma. + + In order for this to work correctly you must provide the correct + values for the camera gain and the camera offset. You can + verify that this is true by loading the histograms and checking + that their shape is approximately that of a Gaussian with 0 + mean and 1 sigma. + """ + def __init__(self, dataSet, parameters=None, analysisName=None): + super().__init__(dataSet, parameters, analysisName) + + if 'highpass_sigma' not in self.parameters: + self.parameters['highpass_sigma'] = 3 + if 'filter_iterations' not in self.parameters: + self.parameters['filter_iterations'] = 5 + + self._cameraGain = self.parameters['camera_gain'] + self._cameraOffset = self.parameters['camera_offset'] + self._filterIterations = self.parameters['filter_iterations'] + self._highPassSigma = self.parameters['highpass_sigma'] + self._highPassFilterSize = int(2 * np.ceil(2 * self._highPassSigma) + 1) + + def _run_analysis(self, fragmentIndex): + warpTask = self.dataSet.load_analysis_task( + self.parameters['warp_task']) + + histogramBins = np.arange(0, 1000, 1) + pixelHistogram = np.zeros( + (self.get_codebook().get_bit_count(), len(histogramBins)-1)) + + # this currently only is to calculate the pixel histograms in order + # to estimate the initial scale factors. This is likely unnecessary + for bi, b in enumerate(self.get_codebook().get_bit_names()): + dataChannel = self.dataSet.get_data_organization()\ + .get_data_channel_for_bit(b) + for i in range(len(self.dataSet.get_z_positions())): + inputImage = warpTask.get_aligned_image( + fragmentIndex, dataChannel, i) + significanceImage = self._preprocess_image(inputImage) + + pixelHistogram[bi, :] += np.histogram( + significanceImage, bins=histogramBins)[0] + + self._save_pixel_histogram(pixelHistogram, fragmentIndex) + + def _preprocess_image(self, inputImage: np.ndarray) -> np.ndarray: + [fg, bg] = imagefilters.high_low_filter(inputImage, + self._highPassFilterSize, + self._highPassSigma, + self._filterIterations) + return imagefilters.est_significance(fg, bg) From 5edfbbbbb648486f69c38c0c168e155331fc64c1 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Fri, 1 May 2020 15:55:27 -0400 Subject: [PATCH 05/12] Initial implementation of 'complete' SNB pipeline, untested. --- merlin/analysis/decode.py | 48 ++++++++++--- merlin/analysis/optimize.py | 65 ++++++++++++++++++ merlin/analysis/preprocess.py | 11 ++- merlin/util/decoding.py | 123 +++++++++++++++++++++++++++++----- 4 files changed, 218 insertions(+), 29 deletions(-) diff --git a/merlin/analysis/decode.py b/merlin/analysis/decode.py index 2ae7e98f..5799cb37 100755 --- a/merlin/analysis/decode.py +++ b/merlin/analysis/decode.py @@ -89,6 +89,13 @@ def get_codebook(self) -> Codebook: self.parameters['preprocess_task']) return preprocessTask.get_codebook() + def _get_decoder(self, codebook, scaleFactors, backgrounds): + decoder = decoding.PixelBasedDecoder(codebook) + return lambda x: \ + decoder.decode_pixels(x, scaleFactors, backgrounds, + lowPassSigma=self.parameters['lowpass_sigma'], + distanceThreshold=self.parameters['distance_threshold']) + def _run_analysis(self, fragmentIndex): """This function decodes the barcodes in a fov and saves them to the barcode database. @@ -102,9 +109,9 @@ def _run_analysis(self, fragmentIndex): lowPassSigma = self.parameters['lowpass_sigma'] codebook = self.get_codebook() - decoder = decoding.PixelBasedDecoder(codebook) scaleFactors = optimizeTask.get_scale_factors() backgrounds = optimizeTask.get_backgrounds() + decodeFn = self._get_decoder(codebook, scaleFactors, backgrounds) chromaticCorrector = optimizeTask.get_chromatic_corrector() zPositionCount = len(self.dataSet.get_z_positions()) @@ -145,10 +152,7 @@ def _run_analysis(self, fragmentIndex): (imageSet.shape[0], imageSet.shape[-2], imageSet.shape[-1])) - di, pm, npt, d = decoder.decode_pixels( - imageSet, scaleFactors, backgrounds, - lowPassSigma=lowPassSigma, - distanceThreshold=self.parameters['distance_threshold']) + di, pm, npt, d = decodeFn(imageSet) normalizedPixelTraces[zIndex, :, :, :] = npt decodedImages[zIndex, :, :] = di @@ -176,18 +180,16 @@ def _run_analysis(self, fragmentIndex): def _process_independent_z_slice( - self, fov: int, zIndex: int, chromaticCorrector, scaleFactors, - backgrounds, preprocessTask, decoder): + self, fov: int, zIndex: int, chromaticCorrector, preprocessTask, + decoderFn): imageSet = preprocessTask.get_processed_image_set( fov, zIndex, chromaticCorrector) imageSet = imageSet.reshape( (imageSet.shape[0], imageSet.shape[-2], imageSet.shape[-1])) - di, pm, npt, d = decoder.decode_pixels( - imageSet, scaleFactors, backgrounds, - lowPassSigma=self.parameters['lowpass_sigma'], - distanceThreshold=self.parameters['distance_threshold']) + di, pm, npt, d = decoderFn(imageSet) + self._extract_and_save_barcodes( decoder, di, pm, npt, d, fov, zIndex) @@ -235,3 +237,27 @@ def _remove_z_duplicate_barcodes(self, bc): self.parameters['z_duplicate_xy_pixel_threshold'], self.dataSet.get_z_positions()) return bc + + +class DecodeSNB(Decode): + + """ + This variant is designed for shot noise based analysis. + """ + + def __init__(self, dataSet: dataset.MERFISHDataSet, + parameters=None, analysisName=None): + super().__init__(dataSet, parameters, analysisName) + + if 'lowpass_sigma' not in parameters: + self.parameters['lowpass_sigma'] = -1 + if 'significance_threshold' not in self.parameters: + self.parameters['significance_threshold'] = 6 + + def _get_decoder(self, codebook, scaleFactors, backgrounds): + decoder = decoding.PixelBasedDecoder(codebook) + return lambda x: \ + decoder.decode_pixels( + x, significanceThreshold=self.parameters['significance_threshold'], + lowPassSigma=self.parameters['lowpass_sigma'], + distanceThreshold=self.parameters['distance_threshold']) diff --git a/merlin/analysis/optimize.py b/merlin/analysis/optimize.py index f182ab2b..98c8ffc8 100755 --- a/merlin/analysis/optimize.py +++ b/merlin/analysis/optimize.py @@ -414,3 +414,68 @@ def get_barcode_count_history(self) -> np.ndarray: self.parameters['previous_iteration'] ).get_barcode_count_history() return np.append(previousHistory, [countsMean], axis=0) + + + +class OptimizeIterationChromaticCorrection(OptimizeIteration): + + """ + This variant only optimizes the chromatic correction. It is + used in the shot noise based analysis pathway where optimizing + scale factors and background offsets isn't relevant. + """ + + def __init__(self, dataSet, parameters=None, analysisName=None): + super().__init__(dataSet, parameters, analysisName) + + if 'significance_threshold' not in self.parameters: + self.parameters['significance_threshold'] = 10 + + def get_backgrounds(self) -> np.ndarray: + return None + + def get_scale_factors(self) -> np.ndarray: + return None + + def _run_analysis(self, fragmentIndex): + if not self.parameters['optimize_chromatic_correction']: + return + + preprocessTask = self.dataSet.load_analysis_task( + self.parameters['preprocess_task']) + codebook = self.get_codebook() + + fovIndex, zIndex = self.parameters['fov_index'][fragmentIndex] + + chromaticTransformations = \ + self._get_previous_chromatic_transformations() + + self.dataSet.save_pickle_analysis_result( + chromaticTransformations, 'previous_chromatic_corrections', + self.analysisName, resultIndex=fragmentIndex) + self.dataSet.save_numpy_analysis_result( + np.array([fovIndex, zIndex]), 'select_frame', self.analysisName, + resultIndex=fragmentIndex) + + chromaticCorrector = aberration.RigidChromaticCorrector( + chromaticTransformations, self.get_reference_color()) + warpedImages = preprocessTask.get_processed_image_set( + fovIndex, zIndex=zIndex, chromaticCorrector=chromaticCorrector) + + decoder = decoding.PixelBasedDecoderSNB(codebook) + areaThreshold = self.parameters['area_threshold'] + decoder.refactorAreaThreshold = areaThreshold + di, pm, npt, d = \ + decoder.decode_pixels(warpedImages, + signicanceThreshold = \ + self.parameters['significance_threshold']) + + # TODO this saves the barcodes under fragment instead of fov + # the barcodedb should be made more general + cropWidth = self.parameters['crop_width'] + self.get_barcode_database().write_barcodes( + pandas.concat([decoder.extract_barcodes_with_index( + i, di, pm, npt, d, fovIndex, cropWidth, + zIndex, minimumArea=areaThreshold) + for i in range(codebook.get_barcode_count())]), + fov=fragmentIndex) diff --git a/merlin/analysis/preprocess.py b/merlin/analysis/preprocess.py index 8c718842..304fed28 100755 --- a/merlin/analysis/preprocess.py +++ b/merlin/analysis/preprocess.py @@ -184,6 +184,8 @@ class EstimatePixelSignificance(Preprocess): def __init__(self, dataSet, parameters=None, analysisName=None): super().__init__(dataSet, parameters, analysisName) + if 'calculate_histograms' not in self.parameters: + self.parameters['calculate_histograms'] = False if 'highpass_sigma' not in self.parameters: self.parameters['highpass_sigma'] = 3 if 'filter_iterations' not in self.parameters: @@ -196,15 +198,18 @@ def __init__(self, dataSet, parameters=None, analysisName=None): self._highPassFilterSize = int(2 * np.ceil(2 * self._highPassSigma) + 1) def _run_analysis(self, fragmentIndex): + if not self.parameters['calculate_histograms']: + return + warpTask = self.dataSet.load_analysis_task( self.parameters['warp_task']) - histogramBins = np.arange(0, 1000, 1) + histogramBins = np.arange(0, 1000, 0.1) pixelHistogram = np.zeros( (self.get_codebook().get_bit_count(), len(histogramBins)-1)) - # this currently only is to calculate the pixel histograms in order - # to estimate the initial scale factors. This is likely unnecessary + # This is only used to calculate the pixel sigma value histograms, + # which might not be that useful. for bi, b in enumerate(self.get_codebook().get_bit_names()): dataChannel = self.dataSet.get_data_organization()\ .get_data_channel_for_bit(b) diff --git a/merlin/util/decoding.py b/merlin/util/decoding.py index 2447917f..76999a3d 100644 --- a/merlin/util/decoding.py +++ b/merlin/util/decoding.py @@ -14,6 +14,19 @@ """ +def filter_images(imageData, lowPassSigma): + if (lowPassSigma > 0): + filteredImages = np.zeros(imageData.shape, dtype=np.float32) + filterSize = int(2 * np.ceil(2 * lowPassSigma) + 1) + for i in range(imageData.shape[0]): + filteredImages[i, :, :] = cv2.GaussianBlur( + imageData[i, :, :], (filterSize, filterSize), lowPassSigma) + else: + filteredImages = imageData.astype(np.float32) + + return filteredImages + + def normalize(x): norm = np.linalg.norm(x) if norm > 0: @@ -90,21 +103,18 @@ def decode_pixels(self, imageData: np.ndarray, if backgrounds is None: backgrounds = self._backgrounds - filteredImages = np.zeros(imageData.shape, dtype=np.float32) - filterSize = int(2 * np.ceil(2 * lowPassSigma) + 1) - for i in range(imageData.shape[0]): - filteredImages[i, :, :] = cv2.GaussianBlur( - imageData[i, :, :], (filterSize, filterSize), lowPassSigma) + scaleFactors = scaleFactors.astype(np.float32) + backgrounds = backgrounds.astype(np.float32) + + filteredImages = filter_images(imageData, lowPassSigma) pixelTraces = np.reshape( - filteredImages, + filteredImages, (filteredImages.shape[0], np.prod(filteredImages.shape[1:]))) - scaledPixelTraces = np.transpose( - np.array([(p-b)/s for p, s, b in zip(pixelTraces, scaleFactors, - backgrounds)])) - pixelMagnitudes = np.array( - [np.linalg.norm(x) for x in scaledPixelTraces], dtype=np.float32) + scaledPixelTraces = np.transpose(pixelTraces).astype(np.float32) + scaledPixelTraces = (scaledPixelTraces - backgrounds)/scaleFactors + pixelMagnitudes = np.linalg.norm(scaledPixelTraces, axis = 1).astype(np.float32) pixelMagnitudes[pixelMagnitudes == 0] = 1 normalizedPixelTraces = scaledPixelTraces/pixelMagnitudes[:, None] @@ -115,10 +125,9 @@ def decode_pixels(self, imageData: np.ndarray, distances, indexes = neighbors.kneighbors( normalizedPixelTraces, return_distance=True) - decodedImage = np.reshape( - np.array([i[0] if d[0] <= distanceThreshold else -1 - for i, d in zip(indexes, distances)], dtype=np.int16), - filteredImages.shape[1:]) + decodedImage = indexes.astype(np.int16) + decodedImage[distances > distanceThreshold] = -1 + decodedImage = np.reshape(decodedImage, filteredImages.shape[1:]) pixelMagnitudes = np.reshape(pixelMagnitudes, filteredImages.shape[1:]) normalizedPixelTraces = np.moveaxis(normalizedPixelTraces, 1, 0) @@ -375,3 +384,87 @@ def _extract_backgrounds( backgroundRefactors = offBitIntensity return backgroundRefactors + + +class PixelBasedDecoderSNB(PixelBasedDecoder): + + """ + This variant is designed for the shot noise based analysis + pathway. In this pathway the pixel values are not AU but + an estimate of the actually significance in units of sigma + of the pixel, so in this decoder there is no scaling or + background offset. + """ + + def decode_pixels(self, imageData: np.ndarray, + distanceThreshold: float=0.5176, + significanceThreshold: float=6, + lowPassSigma: float=-1): + """Assign barcodes to the pixels in the provided image stock. + + Each pixel is assigned to the nearest barcode from the codebook if + the distance between the normalized pixel trace and the barcode is + less than the distance threshold. + + Args: + imageData: input image stack. The first dimension indexes the bit + number and the second and third dimensions contain the + corresponding image. + distanceThreshold: the maximum distance between an assigned pixel + and the nearest barcode. Pixels for which the nearest barcode + is greater than distanceThreshold are left unassigned. + significanceThreshold: the minimum pixel significance to not be + 'noise' and set to zero. + lowPassSigma: standard deviation for the low pass filter that is + applied to the images prior to decoding. Usually this is not + done because the images were not deconvolved in the first + place. + Returns: + Four results are returned as a tuple (decodedImage, pixelMagnitudes, + normalizedPixelTraces, distances). decodedImage is an image + indicating the barcode index assigned to each pixel. Pixels + for which a barcode is not assigned have a value of -1. + pixelMagnitudes is an image where each pixel is the norm of + the pixel trace after scaling by the provided scaleFactors. + normalizedPixelTraces is an image stack containing the + normalized intensities for each pixel. distances is an + image containing the distance for each pixel to the assigned + barcode. + """ + filteredImages = filter_images(imageData, lowPassSigma) + + pixelTraces = np.reshape( + filteredImages, + (filteredImages.shape[0], np.prod(filteredImages.shape[1:]))) + pixelTraces = np.transpose(pixelTraces) + pixelTraces[(pixelTraces < significanceThreshold)] = 0 + + pixelMagnitudes = np.linalg.norm(pixelTraces, axis = 1).astype(np.float32) + pixelMask = (pixelMagnitudes > 0) + + normalizedPixelTraces = pixelTraces[pixelMask,:]/ \ + pixelMagnitudes[pixelMask, None] + + neighbors = NearestNeighbors(n_neighbors=1, algorithm='ball_tree') + neighbors.fit(self._decodingMatrix) + + distances, indexes = neighbors.kneighbors( + normalizedPixelTraces, return_distance=True) + + nptImages = np.zeros(pixelTraces.shape, dtype = np.float32) + nptImages[pixelMask,:] = normalizedPixelTraces + nptImages = np.moveaxis(nptImages, 1, 0) + nptImages = np.reshape(nptImages, filteredImages.shape) + + pixelMask = np.reshape(pixelMask, filteredImages.shape[1:]) + + distanceImage = np.zeros(pixelMask.shape, dtype = np.float32) + distanceImage[pixelMask] = distances[:,0] + + decodedImage = np.zeros(pixelMask.shape, dtype = np.int16) - 1 + decodedImage[pixelMask] = indexes[:,0] + decodedImage[distanceImage > distanceThreshold] = -1 + + pixelMagnitudes = np.reshape(pixelMagnitudes, filteredImages.shape[1:]) + + return decodedImage, pixelMagnitudes, nptImages, distanceImage From dbbcd6a38f4d79af48bd6f465ea09b3fa118157c Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Fri, 1 May 2020 16:23:04 -0400 Subject: [PATCH 06/12] Remove list comprehensions to improve performance. --- merlin/util/decoding.py | 77 +++++++++++++++++++++-------------------- 1 file changed, 40 insertions(+), 37 deletions(-) diff --git a/merlin/util/decoding.py b/merlin/util/decoding.py index 2447917f..6c64e20c 100644 --- a/merlin/util/decoding.py +++ b/merlin/util/decoding.py @@ -25,7 +25,8 @@ def normalize(x): class PixelBasedDecoder(object): def __init__(self, codebook: mcodebook.Codebook, - scaleFactors: np.ndarray=None, backgrounds: np.ndarray=None): + scaleFactors: np.ndarray = None, + backgrounds: np.ndarray = None): self._codebook = codebook self._decodingMatrix = self._calculate_normalized_barcodes() self._barcodeCount = self._decodingMatrix.shape[0] @@ -44,11 +45,11 @@ def __init__(self, codebook: mcodebook.Codebook, self.refactorAreaThreshold = 4 def decode_pixels(self, imageData: np.ndarray, - scaleFactors: np.ndarray=None, - backgrounds: np.ndarray=None, - distanceThreshold: float=0.5176, - magnitudeThreshold: float=1, - lowPassSigma: float=1): + scaleFactors: np.ndarray = None, + backgrounds: np.ndarray = None, + distanceThreshold: float = 0.5176, + magnitudeThreshold: float = 1, + lowPassSigma: float = 1): """Assign barcodes to the pixels in the provided image stock. Each pixel is assigned to the nearest barcode from the codebook if @@ -74,16 +75,16 @@ def decode_pixels(self, imageData: np.ndarray, lowPassSigma: standard deviation for the low pass filter that is applied to the images prior to decoding. Returns: - Four results are returned as a tuple (decodedImage, pixelMagnitudes, - normalizedPixelTraces, distances). decodedImage is an image - indicating the barcode index assigned to each pixel. Pixels - for which a barcode is not assigned have a value of -1. - pixelMagnitudes is an image where each pixel is the norm of - the pixel trace after scaling by the provided scaleFactors. - normalizedPixelTraces is an image stack containing the - normalized intensities for each pixel. distances is an - image containing the distance for each pixel to the assigned - barcode. + Four results are returned as a tuple (decodedImage, + pixelMagnitudes, normalizedPixelTraces, distances). + decodedImage is an image indicating the barcode index assigned + to each pixel. Pixels for which a barcode is not assigned have + a value of -1. pixelMagnitudes is an image where each pixel is + the norm of the pixel trace after scaling by the provided + scaleFactors. normalizedPixelTraces is an image stack + containing the normalized intensities for each pixel. distances + is an image containing the distance for each pixel to the + assigned barcode. """ if scaleFactors is None: scaleFactors = self._scaleFactors @@ -96,15 +97,17 @@ def decode_pixels(self, imageData: np.ndarray, filteredImages[i, :, :] = cv2.GaussianBlur( imageData[i, :, :], (filterSize, filterSize), lowPassSigma) + scaleFactors = scaleFactors.astype(np.float32) + backgrounds = backgrounds.astype(np.float32) + pixelTraces = np.reshape( - filteredImages, + filteredImages, (filteredImages.shape[0], np.prod(filteredImages.shape[1:]))) - scaledPixelTraces = np.transpose( - np.array([(p-b)/s for p, s, b in zip(pixelTraces, scaleFactors, - backgrounds)])) - pixelMagnitudes = np.array( - [np.linalg.norm(x) for x in scaledPixelTraces], dtype=np.float32) + scaledPixelTraces = np.transpose(pixelTraces).astype(np.float32) + scaledPixelTraces = (scaledPixelTraces - backgrounds)/scaleFactors + pixelMagnitudes = \ + np.linalg.norm(scaledPixelTraces, axis=1).astype(np.float32) pixelMagnitudes[pixelMagnitudes == 0] = 1 normalizedPixelTraces = scaledPixelTraces/pixelMagnitudes[:, None] @@ -115,10 +118,9 @@ def decode_pixels(self, imageData: np.ndarray, distances, indexes = neighbors.kneighbors( normalizedPixelTraces, return_distance=True) - decodedImage = np.reshape( - np.array([i[0] if d[0] <= distanceThreshold else -1 - for i, d in zip(indexes, distances)], dtype=np.int16), - filteredImages.shape[1:]) + decodedImage = indexes.astype(np.int16) + decodedImage[distances > distanceThreshold] = -1 + decodedImage = np.reshape(decodedImage, filteredImages.shape[1:]) pixelMagnitudes = np.reshape(pixelMagnitudes, filteredImages.shape[1:]) normalizedPixelTraces = np.moveaxis(normalizedPixelTraces, 1, 0) @@ -133,8 +135,8 @@ def decode_pixels(self, imageData: np.ndarray, def extract_barcodes_with_index( self, barcodeIndex: int, decodedImage: np.ndarray, pixelMagnitudes: np.ndarray, pixelTraces: np.ndarray, - distances: np.ndarray, fov: int, cropWidth: int, zIndex: int = None, - globalAligner=None, minimumArea: int = 0 + distances: np.ndarray, fov: int, cropWidth: int, + zIndex: int = None, globalAligner=None, minimumArea: int = 0 ) -> pandas.DataFrame: """Extract the barcode information from the decoded image for barcodes that were decoded to the specified barcode index. @@ -151,8 +153,8 @@ def extract_barcodes_with_index( distances: an image indicating the distance between the normalized pixel trace and the assigned barcode for each pixel fov: the index of the field of view - cropWidth: the number of pixels around the edge of each image within - which barcodes are excluded from the output list. + cropWidth: the number of pixels around the edge of each image + within which barcodes are excluded from the output list. zIndex: the index of the z position globalAligner: the aligner used for converted to local x,y coordinates to global x,y coordinates @@ -199,7 +201,8 @@ def extract_barcodes_with_index( else: intensityAndCoords = [ - np.array([[y[0], y[1], pixelMagnitudes[y[0], y[1]]] for y in x]) + np.array([[y[0], y[1], + pixelMagnitudes[y[0], y[1]]] for y in x]) for x in allCoords] centroidCoords = np.array( [[(r[:, 0] * (r[:, -1] / r[:, -1].sum())).sum(), @@ -260,16 +263,16 @@ def _calculate_normalized_barcodes( ignoreBlanks: Flag to set if the barcodes corresponding to blanks should be ignored. If True, barcodes corresponding to a name that contains 'Blank' are ignored. - includeErrors: Flag to set if barcodes corresponding to single bit + includeErrors: Flag to set if barcodes corresponding to single bit errors should be added. Returns: A 2d numpy array where each row is a normalized barcode and each column is the corresponding normalized bit value. """ - + barcodeSet = self._codebook.get_barcodes(ignoreBlanks=ignoreBlanks) magnitudes = np.sqrt(np.sum(barcodeSet*barcodeSet, axis=1)) - + if not includeErrors: weightedBarcodes = np.array( [normalize(x) for x, m in zip(barcodeSet, magnitudes)]) @@ -289,7 +292,7 @@ def _calculate_normalized_barcodes( def extract_refactors( self, decodedImage, pixelMagnitudes, normalizedPixelTraces, - extractBackgrounds = False + extractBackgrounds=False ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """Calculate the scale factors that would result in the mean on bit intensity for each bit to be equal. @@ -349,8 +352,8 @@ def _extract_backgrounds( imageSet: the image stack to decode in order to determine the scale factors Returns: - an array of the backgrounds where the i'th entry is the scale factor - for bit i. + an array of the backgrounds where the i'th entry is the scale + factor for bit i. """ sumMinPixelTraces = np.zeros((self._barcodeCount, self._bitCount)) barcodesSeen = np.zeros(self._barcodeCount) From 6208b116e6cf89a36329d8fb4a91ddf42da0863d Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Fri, 1 May 2020 19:03:39 -0400 Subject: [PATCH 07/12] Fix various bugs in the code such that it now works in the MERlin pipeline. Style code. --- merlin/analysis/decode.py | 74 +++++++++++++----------- merlin/analysis/optimize.py | 23 ++++---- merlin/analysis/preprocess.py | 26 ++++----- merlin/util/decoding.py | 105 +++++++++++++++++----------------- 4 files changed, 120 insertions(+), 108 deletions(-) diff --git a/merlin/analysis/decode.py b/merlin/analysis/decode.py index 5799cb37..defa0001 100755 --- a/merlin/analysis/decode.py +++ b/merlin/analysis/decode.py @@ -91,10 +91,13 @@ def get_codebook(self) -> Codebook: def _get_decoder(self, codebook, scaleFactors, backgrounds): decoder = decoding.PixelBasedDecoder(codebook) - return lambda x: \ - decoder.decode_pixels(x, scaleFactors, backgrounds, - lowPassSigma=self.parameters['lowpass_sigma'], - distanceThreshold=self.parameters['distance_threshold']) + return [decoder, + decoder.decode_pixels( + x, scaleFactors, backgrounds, + lowPassSigma= + self.parameters['lowpass_sigma'], + distanceThreshold= + self.parameters['distance_threshold'])] def _run_analysis(self, fragmentIndex): """This function decodes the barcodes in a fov and saves them to the @@ -111,7 +114,9 @@ def _run_analysis(self, fragmentIndex): codebook = self.get_codebook() scaleFactors = optimizeTask.get_scale_factors() backgrounds = optimizeTask.get_backgrounds() - decodeFn = self._get_decoder(codebook, scaleFactors, backgrounds) + [decoder, decodeFn] = self._get_decoder(codebook, + scaleFactors, + backgrounds) chromaticCorrector = optimizeTask.get_chromatic_corrector() zPositionCount = len(self.dataSet.get_z_positions()) @@ -125,9 +130,8 @@ def _run_analysis(self, fragmentIndex): if not decode3d: for zIndex in range(zPositionCount): di, pm, d = self._process_independent_z_slice( - fragmentIndex, zIndex, chromaticCorrector, scaleFactors, - backgrounds, preprocessTask, decoder - ) + fragmentIndex, zIndex, chromaticCorrector, + preprocessTask, decoder, decodeFn) decodedImages[zIndex, :, :] = di magnitudeImages[zIndex, :, :] = pm @@ -178,10 +182,9 @@ def _run_analysis(self, fragmentIndex): bcDB.empty_database(fragmentIndex) bcDB.write_barcodes(bc, fov=fragmentIndex) - def _process_independent_z_slice( self, fov: int, zIndex: int, chromaticCorrector, preprocessTask, - decoderFn): + decoder, decoderFn): imageSet = preprocessTask.get_processed_image_set( fov, zIndex, chromaticCorrector) @@ -199,25 +202,26 @@ def _save_decoded_images(self, fov: int, zPositionCount: int, decodedImages: np.ndarray, magnitudeImages: np.ndarray, distanceImages: np.ndarray) -> None: - imageDescription = self.dataSet.analysis_tiff_description( - zPositionCount, 3) - with self.dataSet.writer_for_analysis_images( - self, 'decoded', fov) as outputTif: - for i in range(zPositionCount): - outputTif.save(decodedImages[i].astype(np.float32), - photometric='MINISBLACK', - metadata=imageDescription) - outputTif.save(magnitudeImages[i].astype(np.float32), - photometric='MINISBLACK', - metadata=imageDescription) - outputTif.save(distanceImages[i].astype(np.float32), - photometric='MINISBLACK', - metadata=imageDescription) + imageDescription = self.dataSet.analysis_tiff_description( + zPositionCount, 3) + with self.dataSet.writer_for_analysis_images( + self, 'decoded', fov) as outputTif: + for i in range(zPositionCount): + outputTif.save(decodedImages[i].astype(np.float32), + photometric='MINISBLACK', + metadata=imageDescription) + outputTif.save(magnitudeImages[i].astype(np.float32), + photometric='MINISBLACK', + metadata=imageDescription) + outputTif.save(distanceImages[i].astype(np.float32), + photometric='MINISBLACK', + metadata=imageDescription) def _extract_and_save_barcodes( - self, decoder: decoding.PixelBasedDecoder, decodedImage: np.ndarray, - pixelMagnitudes: np.ndarray, pixelTraces: np.ndarray, - distances: np.ndarray, fov: int, zIndex: int=None) -> None: + self, decoder: decoding.PixelBasedDecoder, + decodedImage: np.ndarray, pixelMagnitudes: np.ndarray, + pixelTraces: np.ndarray, distances: np.ndarray, + fov: int, zIndex: int = None) -> None: globalTask = self.dataSet.load_analysis_task( self.parameters['global_align_task']) @@ -255,9 +259,13 @@ def __init__(self, dataSet: dataset.MERFISHDataSet, self.parameters['significance_threshold'] = 6 def _get_decoder(self, codebook, scaleFactors, backgrounds): - decoder = decoding.PixelBasedDecoder(codebook) - return lambda x: \ - decoder.decode_pixels( - x, significanceThreshold=self.parameters['significance_threshold'], - lowPassSigma=self.parameters['lowpass_sigma'], - distanceThreshold=self.parameters['distance_threshold']) + decoder = decoding.PixelBasedDecoderSNB(codebook) + return [decoder, + decoder.decode_pixels( + x, + significanceThreshold= + self.parameters['significance_threshold'], + lowPassSigma= + self.parameters['lowpass_sigma'], + distanceThreshold= + self.parameters['distance_threshold'])] diff --git a/merlin/analysis/optimize.py b/merlin/analysis/optimize.py index 98c8ffc8..d94f5192 100755 --- a/merlin/analysis/optimize.py +++ b/merlin/analysis/optimize.py @@ -332,9 +332,10 @@ def get_scale_factors(self) -> np.ndarray: # Don't rescale bits that were never seen refactors[refactors == 0] = 1 - previousFactors = np.array([self.dataSet.load_numpy_analysis_result( - 'previous_scale_factors', self.analysisName, resultIndex=i) - for i in range(self.parameters['fov_per_iteration'])]) + previousFactors = np.array( + [self.dataSet.load_numpy_analysis_result( + 'previous_scale_factors', self.analysisName, resultIndex=i) + for i in range(self.parameters['fov_per_iteration'])]) scaleFactors = np.nanmedian( np.multiply(refactors, previousFactors), axis=0) @@ -364,9 +365,10 @@ def get_backgrounds(self) -> np.ndarray: 'previous_backgrounds', self.analysisName, resultIndex=i) for i in range(self.parameters['fov_per_iteration'])]) - previousFactors = np.array([self.dataSet.load_numpy_analysis_result( - 'previous_scale_factors', self.analysisName, resultIndex=i) - for i in range(self.parameters['fov_per_iteration'])]) + previousFactors = np.array( + [self.dataSet.load_numpy_analysis_result( + 'previous_scale_factors', self.analysisName, resultIndex=i) + for i in range(self.parameters['fov_per_iteration'])]) backgrounds = np.nanmedian(np.add( previousBackgrounds, np.multiply(refactors, previousFactors)), @@ -416,7 +418,6 @@ def get_barcode_count_history(self) -> np.ndarray: return np.append(previousHistory, [countsMean], axis=0) - class OptimizeIterationChromaticCorrection(OptimizeIteration): """ @@ -433,14 +434,14 @@ def __init__(self, dataSet, parameters=None, analysisName=None): def get_backgrounds(self) -> np.ndarray: return None - + def get_scale_factors(self) -> np.ndarray: return None - + def _run_analysis(self, fragmentIndex): if not self.parameters['optimize_chromatic_correction']: return - + preprocessTask = self.dataSet.load_analysis_task( self.parameters['preprocess_task']) codebook = self.get_codebook() @@ -467,7 +468,7 @@ def _run_analysis(self, fragmentIndex): decoder.refactorAreaThreshold = areaThreshold di, pm, npt, d = \ decoder.decode_pixels(warpedImages, - signicanceThreshold = \ + significanceThreshold= self.parameters['significance_threshold']) # TODO this saves the barcodes under fragment instead of fov diff --git a/merlin/analysis/preprocess.py b/merlin/analysis/preprocess.py index 304fed28..ad4e0235 100755 --- a/merlin/analysis/preprocess.py +++ b/merlin/analysis/preprocess.py @@ -14,13 +14,16 @@ class Preprocess(analysistask.ParallelAnalysisTask): """ An abstract class for preparing data for barcode calling. """ - + def __init__(self, dataSet, parameters=None, analysisName=None): super().__init__(dataSet, parameters, analysisName) if 'codebook_index' not in self.parameters: self.parameters['codebook_index'] = 0 + self.warpTask = self.dataSet.load_analysis_task( + self.parameters['warp_task']) + def _image_name(self, fov): destPath = self.dataSet.get_analysis_subdirectory( self.analysisName, subdirectory='preprocessed_images') @@ -64,14 +67,14 @@ def get_processed_image_set( if zIndex is None: return np.array([[self.get_processed_image( fov, self.dataSet.get_data_organization() - .get_data_channel_for_bit(b), zIndex, chromaticCorrector) - for zIndex in range(len(self.dataSet.get_z_positions()))] - for b in self.get_codebook().get_bit_names()]) + .get_data_channel_for_bit(b), zIndex, chromaticCorrector) + for zIndex in range(len(self.dataSet.get_z_positions()))] + for b in self.get_codebook().get_bit_names()]) else: return np.array([self.get_processed_image( fov, self.dataSet.get_data_organization() - .get_data_channel_for_bit(b), zIndex, chromaticCorrector) - for b in self.get_codebook().get_bit_names()]) + .get_data_channel_for_bit(b), zIndex, chromaticCorrector) + for b in self.get_codebook().get_bit_names()]) def get_processed_image( self, fov: int, dataChannel: int, zIndex: int, @@ -103,10 +106,6 @@ def __init__(self, dataSet, parameters=None, analysisName=None): self._deconSigma = self.parameters['decon_sigma'] self._deconIterations = self.parameters['decon_iterations'] - self.warpTask = self.dataSet.load_analysis_task( - self.parameters['warp_task']) - - def _high_pass_filter(self, inputImage: np.ndarray) -> np.ndarray: highPassFilterSize = int(2 * np.ceil(2 * self._highPassSigma) + 1) hpImage = imagefilters.high_pass_filter(inputImage, @@ -175,7 +174,7 @@ class EstimatePixelSignificance(Preprocess): """ Estimates pixel significance in units of sigma. - In order for this to work correctly you must provide the correct + In order for this to work correctly you must provide the correct values for the camera gain and the camera offset. You can verify that this is true by loading the histograms and checking that their shape is approximately that of a Gaussian with 0 @@ -195,12 +194,13 @@ def __init__(self, dataSet, parameters=None, analysisName=None): self._cameraOffset = self.parameters['camera_offset'] self._filterIterations = self.parameters['filter_iterations'] self._highPassSigma = self.parameters['highpass_sigma'] - self._highPassFilterSize = int(2 * np.ceil(2 * self._highPassSigma) + 1) + self._highPassFilterSize =\ + int(2 * np.ceil(2 * self._highPassSigma) + 1) def _run_analysis(self, fragmentIndex): if not self.parameters['calculate_histograms']: return - + warpTask = self.dataSet.load_analysis_task( self.parameters['warp_task']) diff --git a/merlin/util/decoding.py b/merlin/util/decoding.py index 76999a3d..91d7f0be 100644 --- a/merlin/util/decoding.py +++ b/merlin/util/decoding.py @@ -38,7 +38,8 @@ def normalize(x): class PixelBasedDecoder(object): def __init__(self, codebook: mcodebook.Codebook, - scaleFactors: np.ndarray=None, backgrounds: np.ndarray=None): + scaleFactors: np.ndarray = None, + backgrounds: np.ndarray = None): self._codebook = codebook self._decodingMatrix = self._calculate_normalized_barcodes() self._barcodeCount = self._decodingMatrix.shape[0] @@ -57,11 +58,11 @@ def __init__(self, codebook: mcodebook.Codebook, self.refactorAreaThreshold = 4 def decode_pixels(self, imageData: np.ndarray, - scaleFactors: np.ndarray=None, - backgrounds: np.ndarray=None, - distanceThreshold: float=0.5176, - magnitudeThreshold: float=1, - lowPassSigma: float=1): + scaleFactors: np.ndarray = None, + backgrounds: np.ndarray = None, + distanceThreshold: float = 0.5176, + magnitudeThreshold: float = 1, + lowPassSigma: float = 1): """Assign barcodes to the pixels in the provided image stock. Each pixel is assigned to the nearest barcode from the codebook if @@ -87,16 +88,16 @@ def decode_pixels(self, imageData: np.ndarray, lowPassSigma: standard deviation for the low pass filter that is applied to the images prior to decoding. Returns: - Four results are returned as a tuple (decodedImage, pixelMagnitudes, - normalizedPixelTraces, distances). decodedImage is an image - indicating the barcode index assigned to each pixel. Pixels - for which a barcode is not assigned have a value of -1. - pixelMagnitudes is an image where each pixel is the norm of - the pixel trace after scaling by the provided scaleFactors. - normalizedPixelTraces is an image stack containing the - normalized intensities for each pixel. distances is an - image containing the distance for each pixel to the assigned - barcode. + Four results are returned as a tuple (decodedImage, + pixelMagnitudes, normalizedPixelTraces, distances). + decodedImage is an image indicating the barcode index assigned + to each pixel. Pixels for which a barcode is not assigned have + a value of -1. pixelMagnitudes is an image where each pixel is + the norm of the pixel trace after scaling by the provided + scaleFactors. normalizedPixelTraces is an image stack + containing the normalized intensities for each pixel. distances + is an image containing the distance for each pixel to the + assigned barcode. """ if scaleFactors is None: scaleFactors = self._scaleFactors @@ -114,7 +115,8 @@ def decode_pixels(self, imageData: np.ndarray, scaledPixelTraces = np.transpose(pixelTraces).astype(np.float32) scaledPixelTraces = (scaledPixelTraces - backgrounds)/scaleFactors - pixelMagnitudes = np.linalg.norm(scaledPixelTraces, axis = 1).astype(np.float32) + pixelMagnitudes = \ + np.linalg.norm(scaledPixelTraces, axis=1).astype(np.float32) pixelMagnitudes[pixelMagnitudes == 0] = 1 normalizedPixelTraces = scaledPixelTraces/pixelMagnitudes[:, None] @@ -142,8 +144,8 @@ def decode_pixels(self, imageData: np.ndarray, def extract_barcodes_with_index( self, barcodeIndex: int, decodedImage: np.ndarray, pixelMagnitudes: np.ndarray, pixelTraces: np.ndarray, - distances: np.ndarray, fov: int, cropWidth: int, zIndex: int = None, - globalAligner=None, minimumArea: int = 0 + distances: np.ndarray, fov: int, cropWidth: int, + zIndex: int = None, globalAligner=None, minimumArea: int = 0 ) -> pandas.DataFrame: """Extract the barcode information from the decoded image for barcodes that were decoded to the specified barcode index. @@ -160,8 +162,8 @@ def extract_barcodes_with_index( distances: an image indicating the distance between the normalized pixel trace and the assigned barcode for each pixel fov: the index of the field of view - cropWidth: the number of pixels around the edge of each image within - which barcodes are excluded from the output list. + cropWidth: the number of pixels around the edge of each image + within which barcodes are excluded from the output list. zIndex: the index of the z position globalAligner: the aligner used for converted to local x,y coordinates to global x,y coordinates @@ -208,7 +210,8 @@ def extract_barcodes_with_index( else: intensityAndCoords = [ - np.array([[y[0], y[1], pixelMagnitudes[y[0], y[1]]] for y in x]) + np.array([[y[0], y[1], + pixelMagnitudes[y[0], y[1]]] for y in x]) for x in allCoords] centroidCoords = np.array( [[(r[:, 0] * (r[:, -1] / r[:, -1].sum())).sum(), @@ -269,16 +272,16 @@ def _calculate_normalized_barcodes( ignoreBlanks: Flag to set if the barcodes corresponding to blanks should be ignored. If True, barcodes corresponding to a name that contains 'Blank' are ignored. - includeErrors: Flag to set if barcodes corresponding to single bit + includeErrors: Flag to set if barcodes corresponding to single bit errors should be added. Returns: A 2d numpy array where each row is a normalized barcode and each column is the corresponding normalized bit value. """ - + barcodeSet = self._codebook.get_barcodes(ignoreBlanks=ignoreBlanks) magnitudes = np.sqrt(np.sum(barcodeSet*barcodeSet, axis=1)) - + if not includeErrors: weightedBarcodes = np.array( [normalize(x) for x, m in zip(barcodeSet, magnitudes)]) @@ -298,7 +301,7 @@ def _calculate_normalized_barcodes( def extract_refactors( self, decodedImage, pixelMagnitudes, normalizedPixelTraces, - extractBackgrounds = False + extractBackgrounds=False ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """Calculate the scale factors that would result in the mean on bit intensity for each bit to be equal. @@ -358,8 +361,8 @@ def _extract_backgrounds( imageSet: the image stack to decode in order to determine the scale factors Returns: - an array of the backgrounds where the i'th entry is the scale factor - for bit i. + an array of the backgrounds where the i'th entry is the scale + factor for bit i. """ sumMinPixelTraces = np.zeros((self._barcodeCount, self._bitCount)) barcodesSeen = np.zeros(self._barcodeCount) @@ -397,9 +400,9 @@ class PixelBasedDecoderSNB(PixelBasedDecoder): """ def decode_pixels(self, imageData: np.ndarray, - distanceThreshold: float=0.5176, - significanceThreshold: float=6, - lowPassSigma: float=-1): + distanceThreshold: float = 0.5176, + significanceThreshold: float = 6, + lowPassSigma: float = -1): """Assign barcodes to the pixels in the provided image stock. Each pixel is assigned to the nearest barcode from the codebook if @@ -420,30 +423,30 @@ def decode_pixels(self, imageData: np.ndarray, done because the images were not deconvolved in the first place. Returns: - Four results are returned as a tuple (decodedImage, pixelMagnitudes, - normalizedPixelTraces, distances). decodedImage is an image - indicating the barcode index assigned to each pixel. Pixels - for which a barcode is not assigned have a value of -1. - pixelMagnitudes is an image where each pixel is the norm of - the pixel trace after scaling by the provided scaleFactors. - normalizedPixelTraces is an image stack containing the - normalized intensities for each pixel. distances is an - image containing the distance for each pixel to the assigned - barcode. + Four results are returned as a tuple (decodedImage, + pixelMagnitudes, normalizedPixelTraces, distances). + decodedImage is an image indicating the barcode index assigned + to each pixel. Pixels for which a barcode is not assigned have + a value of -1. pixelMagnitudes is an image where each pixel is + the norm of the pixel trace after scaling by the provided + scaleFactors. normalizedPixelTraces is an image stack + containing the normalized intensities for each pixel. distances + is an image containing the distance for each pixel to the """ filteredImages = filter_images(imageData, lowPassSigma) pixelTraces = np.reshape( - filteredImages, + filteredImages, (filteredImages.shape[0], np.prod(filteredImages.shape[1:]))) pixelTraces = np.transpose(pixelTraces) pixelTraces[(pixelTraces < significanceThreshold)] = 0 - pixelMagnitudes = np.linalg.norm(pixelTraces, axis = 1).astype(np.float32) + pixelMagnitudes = \ + np.linalg.norm(pixelTraces, axis=1).astype(np.float32) pixelMask = (pixelMagnitudes > 0) - normalizedPixelTraces = pixelTraces[pixelMask,:]/ \ - pixelMagnitudes[pixelMask, None] + normalizedPixelTraces = pixelTraces[pixelMask, :] / \ + pixelMagnitudes[pixelMask, None] neighbors = NearestNeighbors(n_neighbors=1, algorithm='ball_tree') neighbors.fit(self._decodingMatrix) @@ -451,18 +454,18 @@ def decode_pixels(self, imageData: np.ndarray, distances, indexes = neighbors.kneighbors( normalizedPixelTraces, return_distance=True) - nptImages = np.zeros(pixelTraces.shape, dtype = np.float32) - nptImages[pixelMask,:] = normalizedPixelTraces + nptImages = np.zeros(pixelTraces.shape, dtype=np.float32) + nptImages[pixelMask, :] = normalizedPixelTraces nptImages = np.moveaxis(nptImages, 1, 0) nptImages = np.reshape(nptImages, filteredImages.shape) pixelMask = np.reshape(pixelMask, filteredImages.shape[1:]) - distanceImage = np.zeros(pixelMask.shape, dtype = np.float32) - distanceImage[pixelMask] = distances[:,0] + distanceImage = np.zeros(pixelMask.shape, dtype=np.float32) + distanceImage[pixelMask] = distances[:, 0] - decodedImage = np.zeros(pixelMask.shape, dtype = np.int16) - 1 - decodedImage[pixelMask] = indexes[:,0] + decodedImage = np.zeros(pixelMask.shape, dtype=np.int16) - 1 + decodedImage[pixelMask] = indexes[:, 0] decodedImage[distanceImage > distanceThreshold] = -1 pixelMagnitudes = np.reshape(pixelMagnitudes, filteredImages.shape[1:]) From 874bea6f4474473cdb7e2c4c194c3e78ba0ad274 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Fri, 1 May 2020 19:21:57 -0400 Subject: [PATCH 08/12] Fix bug introduced by trying to conform with a pycodestyle suggestion. --- merlin/analysis/decode.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/merlin/analysis/decode.py b/merlin/analysis/decode.py index defa0001..f474a48a 100755 --- a/merlin/analysis/decode.py +++ b/merlin/analysis/decode.py @@ -92,7 +92,7 @@ def get_codebook(self) -> Codebook: def _get_decoder(self, codebook, scaleFactors, backgrounds): decoder = decoding.PixelBasedDecoder(codebook) return [decoder, - decoder.decode_pixels( + lambda x: decoder.decode_pixels( x, scaleFactors, backgrounds, lowPassSigma= self.parameters['lowpass_sigma'], @@ -261,7 +261,7 @@ def __init__(self, dataSet: dataset.MERFISHDataSet, def _get_decoder(self, codebook, scaleFactors, backgrounds): decoder = decoding.PixelBasedDecoderSNB(codebook) return [decoder, - decoder.decode_pixels( + lambda x:decoder.decode_pixels( x, significanceThreshold= self.parameters['significance_threshold'], From 3a674f65ce4b5e361370174518efff84bcc689fc Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Sat, 2 May 2020 09:25:04 -0400 Subject: [PATCH 09/12] Actually apply camera parameters to image. --- merlin/analysis/preprocess.py | 1 + 1 file changed, 1 insertion(+) diff --git a/merlin/analysis/preprocess.py b/merlin/analysis/preprocess.py index ad4e0235..fb60a3ce 100755 --- a/merlin/analysis/preprocess.py +++ b/merlin/analysis/preprocess.py @@ -224,6 +224,7 @@ def _run_analysis(self, fragmentIndex): self._save_pixel_histogram(pixelHistogram, fragmentIndex) def _preprocess_image(self, inputImage: np.ndarray) -> np.ndarray: + inputImage = (inputImage - self._cameraOffset)/self._cameraGain [fg, bg] = imagefilters.high_low_filter(inputImage, self._highPassFilterSize, self._highPassSigma, From 2f5a976f31c320ffbad5aa8f820b5af946abe2d4 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Mon, 4 May 2020 12:19:58 -0400 Subject: [PATCH 10/12] Fix image sigma calculation. This now returns the expected answer on simulated images, a Gaussian with zero offset and sigma 1.0. --- merlin/util/imagefilters.py | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/merlin/util/imagefilters.py b/merlin/util/imagefilters.py index 9e05e789..3492b704 100755 --- a/merlin/util/imagefilters.py +++ b/merlin/util/imagefilters.py @@ -21,15 +21,7 @@ def est_significance(foreground: np.ndarray, # a negative number. # background = np.clip(background, 1.0, None) - snbfIm = foreground/np.sqrt(background) - - # Why subtract 1.0 in order to make offset = 0.0? This is - # emperical, not entirely sure I understand why the offset - # is 1.0. - # - snbfIm -= 1.0 - - return snbfIm + return foreground/np.sqrt(background) def high_low_filter(image: np.ndarray, From 7b78c42ae5a52ba2e541f5aff97043da14866878 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Mon, 4 May 2020 12:21:21 -0400 Subject: [PATCH 11/12] Change histogram range. Convert image type to np.float32 to reduce possible rounding errors. --- merlin/analysis/preprocess.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/merlin/analysis/preprocess.py b/merlin/analysis/preprocess.py index fb60a3ce..6c91da4c 100755 --- a/merlin/analysis/preprocess.py +++ b/merlin/analysis/preprocess.py @@ -204,7 +204,7 @@ def _run_analysis(self, fragmentIndex): warpTask = self.dataSet.load_analysis_task( self.parameters['warp_task']) - histogramBins = np.arange(0, 1000, 0.1) + histogramBins = np.arange(-10, 100, 1.0) pixelHistogram = np.zeros( (self.get_codebook().get_bit_count(), len(histogramBins)-1)) @@ -224,6 +224,7 @@ def _run_analysis(self, fragmentIndex): self._save_pixel_histogram(pixelHistogram, fragmentIndex) def _preprocess_image(self, inputImage: np.ndarray) -> np.ndarray: + inputImage = inputImage.astype(np.float32) inputImage = (inputImage - self._cameraOffset)/self._cameraGain [fg, bg] = imagefilters.high_low_filter(inputImage, self._highPassFilterSize, From 99d5bbd4055ffc107af84ff68aa5298a2dbae07a Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Mon, 4 May 2020 13:36:46 -0400 Subject: [PATCH 12/12] Add comment about recommended threshold settings. --- merlin/analysis/decode.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/merlin/analysis/decode.py b/merlin/analysis/decode.py index f474a48a..7e8b57ff 100755 --- a/merlin/analysis/decode.py +++ b/merlin/analysis/decode.py @@ -247,6 +247,10 @@ class DecodeSNB(Decode): """ This variant is designed for shot noise based analysis. + + A significance_threshold value of at least 5.0 is + recommended. A value of 3.0 or 4.0 will return lots of + false positives. """ def __init__(self, dataSet: dataset.MERFISHDataSet,