diff --git a/merlin/analysis/decode.py b/merlin/analysis/decode.py index 2ae7e98f..7e8b57ff 100755 --- a/merlin/analysis/decode.py +++ b/merlin/analysis/decode.py @@ -89,6 +89,16 @@ 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 [decoder, + 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 +112,11 @@ 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() + [decoder, decodeFn] = self._get_decoder(codebook, + scaleFactors, + backgrounds) chromaticCorrector = optimizeTask.get_chromatic_corrector() zPositionCount = len(self.dataSet.get_z_positions()) @@ -118,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 @@ -145,10 +156,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 @@ -174,20 +182,17 @@ 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, scaleFactors, - backgrounds, preprocessTask, decoder): + self, fov: int, zIndex: int, chromaticCorrector, preprocessTask, + decoder, 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) @@ -197,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']) @@ -235,3 +241,35 @@ 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. + + 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, + 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.PixelBasedDecoderSNB(codebook) + return [decoder, + 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..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)), @@ -414,3 +416,67 @@ 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, + significanceThreshold= + 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 32a904f9..6c91da4c 100755 --- a/merlin/analysis/preprocess.py +++ b/merlin/analysis/preprocess.py @@ -15,52 +15,24 @@ class Preprocess(analysistask.ParallelAnalysisTask): An abstract class for preparing data for barcode calling. """ - 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 _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 _save_pixel_histogram(self, histogram, fov): + self.dataSet.save_numpy_analysis_result( + histogram, 'pixel_histogram', self.analysisName, fov, 'histograms') + def fragment_count(self): return len(self.dataSet.get_fovs()) @@ -76,6 +48,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 @@ -83,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, @@ -100,6 +84,28 @@ 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'] + 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 +168,66 @@ 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 '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: + 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): + if not self.parameters['calculate_histograms']: + return + + warpTask = self.dataSet.load_analysis_task( + self.parameters['warp_task']) + + histogramBins = np.arange(-10, 100, 1.0) + pixelHistogram = np.zeros( + (self.get_codebook().get_bit_count(), len(histogramBins)-1)) + + # 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) + 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: + inputImage = inputImage.astype(np.float32) + inputImage = (inputImage - self._cameraOffset)/self._cameraGain + [fg, bg] = imagefilters.high_low_filter(inputImage, + self._highPassFilterSize, + self._highPassSigma, + self._filterIterations) + return imagefilters.est_significance(fg, bg) diff --git a/merlin/util/decoding.py b/merlin/util/decoding.py index 2447917f..719d4cc3 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: @@ -25,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] @@ -44,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 @@ -74,37 +88,38 @@ 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 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) + + 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 +130,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 +147,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 +165,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 +213,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 +275,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 +304,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 +364,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) @@ -375,3 +390,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 + """ + 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 diff --git a/merlin/util/imagefilters.py b/merlin/util/imagefilters.py index 15e62185..3492b704 100755 --- a/merlin/util/imagefilters.py +++ b/merlin/util/imagefilters.py @@ -6,6 +6,51 @@ """ +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). + """ + # 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) + return foreground/np.sqrt(background) + + +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 +68,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 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",