Merge pull request #13 from zengjiexia/SR-analysis-development

Super resolution image analysis development

README.md

@@ -6,9 +6,9 @@ A platform for the rapid analysis of both diffraction-limited and single-molecul
 
 Applications
 ------------
-- Diffraction-Limited Analysis (DFLSP)
+- Diffraction-Limited Analysis
 
-	This program is intended for rapid analysis of diffraction-limited images, to distinguish the protein aggregates (or other targets) captured by the SiMPull surface. The analysis process can also be used for other diffraction-limited images.
+	This program is intended for the rapid analysis of diffraction-limited images generated under fluorescence microscope. It is capable of distinguishing and characterising the protein aggregates (or other targets) captured on surface. 
 
 	Two methods are available:
 		
@@ -17,12 +17,20 @@ Applications
 	2. Trevor - A pure python alternative for the ComDet. (Written by Yunzhao Wu)
 
 	An advanced thresholding method, Orthogonal Analysis, is provided. It can help you to distinguish the actual particles detected from the backgroud noises by setting a threshold based on the 'intensity per area' distribution of the spots.
+
+
+- Super-resolution Image Analysis
+
+	Super-resolution images are stack of images taken when the fluorescent indicator being constantly bound and dissociated from the targets captured on the surface. This feature provides two image reconstruction methods - [GDSC SMLM 1](https://gdsc-smlm.readthedocs.io/en/latest/) and [ThunderSTORM](https://zitmen.github.io/thunderstorm/) to track the fluorescent spot through the stack and reconstruct into a image with higher resolution than diffraction limit.
+
+	After image reconstruction, drift correction and particle clustering methods are also provided in the program, allowing you to further characterise the targets detected.
 
 
 - Liposome Assay Analysis
 
-	This workflow is the improved version of [Calcium Influx Assay](https://github.com/zengjiexia/CalciumInfluxAssay). It is an automated analysis program for [Ultrasensitive Measurement of Ca2+ Influx into Lipid Vesicles Induced by Protein Aggregates](https://doi.org/10.1002/anie.201700966) developed in the Klenerman Group. 
-
+	This workflow is the improved version of [Calcium Influx Assay](https://github.com/zengjiexia/CalciumInfluxAssay). It is an automated analysis program for [Ultrasensitive Measurement of Ca2+ Influx into Lipid Vesicles Induced by Protein Aggregates](https://doi.org/10.1002/anie.201700966) developed at the Klenerman Group. 
+
+
 
 Requirements
 ------------
@@ -31,6 +39,7 @@ Requirements
 	- os
 	- sys
 	- re
+	- datatime
 	- [opencv - 4.1.2.30](https://pypi.org/project/opencv-contrib-python/)
 	- [tqdm - 4.60.0](https://pypi.org/project/tqdm/)
 	- [astropy - 4.0.2](https://www.astropy.org/)
@@ -39,8 +48,9 @@ Requirements
 	- [numpy](https://numpy.org/)
 	- [math](https://docs.python.org/3/library/math.html)
 	- [tifffile](https://pypi.org/project/tifffile/)
-	- [pandas - 0.25.3](https://pandas.pydata.org/)
+	- [pandas - 1.4.1](https://pandas.pydata.org/)
 	- [scikit-image - 0.17.2](https://scikit-image.org/)
+	- [scikit-learn - 0.24.2](https://scikit-learn.org/)
 	- [PySide 6 - 6.0.2](https://pypi.org/project/PySide6/)
 	- [pyqtgraph](https://github.com/pyqtgraph/pyqtgraph)(Included locally)
 	- [pathos](https://pypi.org/project/pathos/)
@@ -50,6 +60,8 @@ Requirements
 
 - [Fiji(is just imagej)](https://imagej.net/Fiji) 
 	- [ComDet](https://github.com/ekatrukha/ComDet)
+	- [GDSC SMLM 1](https://gdsc-smlm.readthedocs.io/en/latest/)
+	- [ThunderSTORM](https://zitmen.github.io/thunderstorm/)
 
 
 Installation

SR_toolkit-master/BinLocalizations.m

@@ -0,0 +1,31 @@
+function [ im ] = BinLocalizations( positions, imsize, zoomfactor )
+
+% Offset and rescale positions
+positions = single(positions);
+positions = positions*zoomfactor;
+%positions = positions + 0.5;
+
+% Filter localizations outside the FOV
+binimsize = ceil(imsize*zoomfactor);
+keep = positions(:,1)>=0;
+keep = keep & positions(:,2)>=0;
+keep = keep & positions(:,1)<=binimsize;
+keep = keep & positions(:,2)<=binimsize;
+positions = positions(keep,:);
+
+% Bin localizations
+im = zeros(binimsize,binimsize);
+for i = 1:size(positions,1)
+    x = round(positions(i,1));
+    y = round(positions(i,2));
+    if x==0
+        x=x+1;
+    end
+    if y==0
+        y=y+1;
+    end
+    im(y,x) = im(y,x)+1;
+end
+
+end
+

SR_toolkit-master/CrossCorrelation.m

@@ -0,0 +1,6 @@
+function [ corr ] = CrossCorrelation( image1, image2 )
+
+    corr = fftshift( real( ifft2(fft2(image1).*conj(fft2(image2))) ) );
+
+end
+

SR_toolkit-master/DCC.m

@@ -0,0 +1,61 @@
+% -------------------------------------------------------------------------
+% Matlab version of direct cross-correlation drift correction method
+% Input:    coords:             localization coordinates [x y frame], 
+%           segpara:            segmentation parameters (time wimdow, frame)
+%           imsize:             image size (pixel)
+%           pixelsize:          camera pixel size (nm)
+%           binsize:            spatial bin pixel size (nm)
+% Output:   coordscorr:         localization coordinates after correction [xc yc] 
+%           finaldrift:         drift curve 
+% By Yina Wang @ Hust 2013.10.18
+% -------------------------------------------------------------------------
+
+function [coordscorr, finaldrift] = DCC(coords, segpara, imsize, pixelsize, binsize )
+
+zoomfactor = pixelsize / binsize;
+ntotalframe = max(coords(:,3));
+
+coordscorr = zeros(size(coords,1),size(coords,2));
+finaldrift = zeros(ntotalframe,2);
+
+%% bin localizations, frame shifts calculation
+nbinframe = floor(ntotalframe / segpara);
+
+imshift = zeros(nbinframe-1,1);
+fitresult = zeros(5,nbinframe-1);
+
+i=1;
+index = coords(:,3)>(i-1)*segpara+1 & coords(:,3)<=i*segpara;
+imorigin = BinLocalizations(coords(index,1:2), imsize, zoomfactor);
+autocorr = CrossCorrelation(imorigin,imorigin);
+[yorigin, xorigin] = GaussianFit(autocorr);
+
+for i = 2:nbinframe
+    index = coords(:,3)>(i-1)*segpara+1 & coords(:,3)<=i*segpara;
+    imbim = BinLocalizations(coords(index,1:2), imsize, zoomfactor);
+
+    corr = CrossCorrelation(imorigin,imbim);
+    [y,x] = GaussianFit(corr);
+    imshift(i-1,1) = xorigin-x;
+    imshift(i-1,2) = yorigin-y;
+end
+
+%% drift interpolation
+indexinterp=zeros(nbinframe+2,1);
+indexinterp(1)=1;
+indexinterp(nbinframe+2)=ntotalframe;
+indexinterp(2:nbinframe+1)=round(segpara/2):segpara:segpara*nbinframe-1;
+
+finaldrift(:,1) = interp1(indexinterp,[0 0 imshift(:,1)' imshift(nbinframe-1,1)],1:ntotalframe,'spline');
+finaldrift(:,2) = interp1(indexinterp,[0 0 imshift(:,2)' imshift(nbinframe-1,2)],1:ntotalframe,'spline');
+finaldrift = finaldrift./zoomfactor;
+
+for i = 1:ntotalframe
+    index = find(coords(:,3)==i);
+    coordscorr(index,1) = coords(index,1)-finaldrift(i,1);
+    coordscorr(index,2) = coords(index,2)-finaldrift(i,2);
+end
+coordscorr(:,3) = coords(:,3);
+
+end
+

SR_toolkit-master/FitInitialGuess.m

@@ -0,0 +1,22 @@
+function [ guess ] = FitInitialGuess( corrfunc, peakfindsize )
+
+    n = size(corrfunc,1)/2;
+
+    [y, x] = ZeroPaddingFFT(corrfunc,peakfindsize,10);
+    guess(4)=int16(ceil(y));
+    guess(3)=int16(ceil(x));
+
+    guess(1)=corrfunc(guess(4),guess(3));
+    guess(5) = mean(mean(corrfunc)); % background
+    cropsize = 64;
+    corr = corrfunc(guess(4)-cropsize/2+1:guess(4)+cropsize/2,guess(3)-cropsize/2+1:guess(3)+cropsize/2);
+    [Wy, Wx] = find( ismember( abs((corr./double(guess(1)) - exp(-1))), min(min(abs(corr./double(guess(1)) - exp(-1)))) ) );
+    Wx = mod(Wx,size(corr,2));
+    guess(2) = mean(( (Wx - cropsize/2 ).^2  + (Wy - cropsize/2 ).^2   ).^(1/2));
+    if(guess(2)<1 || guess(2)>10)
+        guess(2) = 5;
+    end
+
+    guess=double(guess);
+end
+

SR_toolkit-master/GaussianFit.m

@@ -0,0 +1,21 @@
+function [ y, x ] = GaussianFit( corrfunc )
+
+    [ guess ] = FitInitialGuess( corrfunc, 48 );
+
+    radiu = round( 2.5 * guess(2) );
+    originy = guess(4);
+    originx = guess(3);
+    fitdata = corrfunc(originy-radiu:originy+radiu,originx-radiu:originx+radiu);
+    [X,Y] = meshgrid(1:1:size(fitdata,2),1:1:size(fitdata,1));
+    grid = [X Y];
+
+    xoffset = guess(3) - radiu - 1;
+    yoffset = guess(4) - radiu - 1;
+    guess(3) = radiu + 1;
+    guess(4) = radiu + 1;
+    [opt,norm] = lsqcurvefit(@gaussmodel2dwxy, guess, grid, fitdata);
+
+    x = opt(3) + xoffset;
+    y = opt(4) + yoffset;
+end
+

SR_toolkit-master/MCC.m

@@ -0,0 +1,84 @@
+
+% -------------------------------------------------------------------------
+% Matlab version of Mean cross-correction method of OE2012 paper
+% Input:    coords:             localization coordinates [x y frame], 
+%           segpara:            segmentation parameters (time wimdow, frame)
+%           imsize:             image size (pixel)
+%           pixelsize:          camera pixel size (nm)
+%           binsize:            spatial bin pixel size (nm)
+% Output:   coordscorr:         localization coordinates after correction [xc yc] 
+%           finaldrift:         drift curve (save b matrix for other analysis)
+% By Yina Wang @ Hust 2013.10.21
+% -------------------------------------------------------------------------
+
+function [coordscorr, finaldrift, bsave] = MCC(coords, segpara, imsize, pixelsize, binsize )
+
+zoomfactor = pixelsize / binsize;
+ntotalframe = max(coords(:,3));
+
+coordscorr = zeros(size(coords,1),size(coords,2));
+finaldrift = zeros(ntotalframe,2);
+
+%% bin localizations, frame to frame drift calculation
+nbinframe = floor(ntotalframe / segpara);
+
+imshift = zeros(nbinframe * (nbinframe - 1)/2,2); % shifts between any two bin steps in x and y direction
+drift=zeros(nbinframe,2);
+
+flag=1;
+for i = 1:nbinframe-1
+    index = coords(:,3)>(i-1)*segpara+1 & coords(:,3)<=i*segpara;
+    imorigin = BinLocalizations(coords(index,1:2), imsize, zoomfactor);
+    autocorr = CrossCorrelation(imorigin,imorigin);
+    [yorigin, xorigin] = GaussianFit(autocorr);
+
+    for j = i+1:nbinframe
+        index = coords(:,3)>(j-1)*segpara+1 & coords(:,3)<=j*segpara;
+        imbin = BinLocalizations(coords(index,1:2), imsize, zoomfactor);
+        corr = CrossCorrelation(imorigin,imbin);
+        [y,x] = GaussianFit(corr);
+
+        imshift(flag,1) = xorigin-x;
+        imshift(flag,2) = yorigin-y;
+        drift(i,1)=drift(i,1)-imshift(flag,1);
+        drift(i,2)=drift(i,2)-imshift(flag,2);
+        flag=flag+1;
+    end
+
+    clear index
+    for j = 1:i-1
+        index=(nbinframe-1)*(j-1)-(j-2)*(j-1)/2-j+i;
+        drift(i,1)=drift(i,1)+imshift(index,1);
+        drift(i,2)=drift(i,2)+imshift(index,2);
+    end
+end
+for j = 1:nbinframe-1
+    index=(nbinframe-1)*(j-1)-(j-2)*(j-1)/2-j+nbinframe;
+    drift(nbinframe,1)=drift(nbinframe,1)+imshift(index,1);
+    drift(nbinframe,2)=drift(nbinframe,2)+imshift(index,2);
+end
+drift=drift./nbinframe;
+
+%% save results
+bsave=imshift;
+
+%% drift interpolation
+indexinterp=zeros(nbinframe+2,1);
+indexinterp(1)=1;
+indexinterp(nbinframe+2)=ntotalframe;
+indexinterp(2:nbinframe+1)=round(segpara/2):segpara:segpara*nbinframe-1;
+
+drift(:,1)=drift(:,1)-drift(1,1);
+drift(:,2)=drift(:,2)-drift(1,2);
+finaldrift(:,1) = interp1(indexinterp,[0 drift(:,1)' drift(nbinframe,1)],1:ntotalframe,'spline');
+finaldrift(:,2) = interp1(indexinterp,[0 drift(:,2)' drift(nbinframe,2)],1:ntotalframe,'spline');
+finaldrift = finaldrift./zoomfactor;
+
+for i = 1:ntotalframe
+    index = find(coords(:,3)==i);
+    coordscorr(index,1) = coords(index,1)-finaldrift(i,1);
+    coordscorr(index,2) = coords(index,2)-finaldrift(i,2);
+end
+coordscorr(:,3) = coords(:,3);
+
+end