%% Generate 2D gaussian
clear
close all

% We can simulate a gaussian funciton using the fspecial command that we
% saw before to make Gaussian filters
gauss2d = 20000*fspecial('Gaussian',200,20);
windowSize = length(gauss2d);
windowLim = round(windowSize/2);

figure
surf(gauss2d)


figure
surf(gauss2d,'EdgeColor','none')

% Now lets put the Gaussian somewhere in a 5000x5000 image.  We start by
% using randi to generate random integers
nPixels = 5000;
xVal = randi(nPixels,1);
yVal = randi(nPixels,1);

% We'll first make a big image of zeros
img = zeros(nPixels,nPixels);

% Then we'll add the Gaussian filter at the random spot
img(yVal-windowLim:yVal+windowLim-1,xVal-windowLim:xVal+windowLim-1) = gauss2d;

figure
imshow(img)

% Now lets truncate the image size, so the value will be some intermediate
% pixel range.  We can use imresize.  It asks for the image, and then the
% scaling factor
imgTrunc = imresize(img,.1);

figure
imshow(imgTrunc,[])

% Nowe we can add some background noise
imgTruncNoise = mat2gray(imgTrunc + 2 * rand(500));

figure
imshow(imgTruncNoise,[])
% This is our final simulated image.  It contains 1 gaussian fluorophore
% that we can try to fit

clear
close all

gauss2d = 20000*fspecial('Gaussian',200,20);
windowSize = length(gauss2d);
windowLim = round(windowSize/2);

nPixels = 5000;
xVal = randi(nPixels,1);
yVal = randi(nPixels,1);

img = zeros(nPixels,nPixels);
img(yVal-windowLim:yVal+windowLim-1,xVal-windowLim:xVal+windowLim-1) = gauss2d;
imgTrunc = imresize(img,.1);
imgTruncNoise = mat2gray(imgTrunc + 2 * rand(500));

figure
imshow(imgTruncNoise,[])


%% Fitting 2d gaussian

% If we know that our PSF has a width of 2 pixels (assuming the Nyquist
% sampling rate, we can make a Gaussian filter to increase the signal.
gFilt = fspecial('Gaussian',7,2);

% We can convolve the image with our PSF gaussian estimate
imgFilt = mat2gray(conv2(imgTruncNoise,gFilt,'same'));

figure
imshow([imgTruncNoise,imgFilt],[])

% Now we can use an im2bw and regionprops exactly like we've done before
imgMask = im2bw(imgFilt,.3);

figure
imshow([imgFilt imgMask])

% The bounding box is the smallest rectangle that completely encoses the
% region selected by regionprops
stats = regionprops(imgMask,'BoundingBox');
boundbox = stats.BoundingBox

xPixels = floor(boundbox(1))-2:floor(boundbox(1))+boundbox(3)+2;
yPixels = floor(boundbox(2))-2:floor(boundbox(2))+boundbox(4)+2;
subImg = imgTruncNoise(round(yPixels),round(xPixels));

figure
imshow(subImg,'InitialMagnification','Fit')

% Now we just have to fit the 2d gaussian.  First lets define our function
g2DFun = @(a1,sigmax,sigmay,x0,y0,x,y)a1*exp(-((x-x0).^2/(2*sigmax^2) + (y-y0).^2/(2*sigmay^2)));
g2DFitType = fittype(g2DFun,'independent',{'x','y'},'dependent',{'z'});

%Initial guesses
a1 = max(subImg(:));
sigmax = 2;
sigmay = 2;
x0 = length(xPixels)/2;
y0 = length(yPixels)/2;

% Define grid for Matlab to fit.  This just defines the x,y, and z axes to
% fit against.
[X,Y] = meshgrid(1:length(yPixels),1:length(xPixels));
X = X(:);
Y = Y(:);
Z = subImg(:);
fittedData = fit([X,Y],double(Z),g2DFitType,'StartPoint',[a1, sigmax, sigmay, x0, y0])

% The fitetd data variable is like a structure; we can pull out different
% components using the fittedData.Component syntax
x0Fit = fittedData.x0 + xPixels(1)-1;
y0Fit = fittedData.y0 + yPixels(1)-1;
xVal/10 - x0Fit
yVal/10 - y0Fit
figure
s1 = plot(fittedData);
alpha(s1,.4);
hold on
s2 = surf(subImg);
alpha(s2,.6);
hold off


% Lets look at residuals
fitDataTrunc = g2DFun(fittedData.a1, fittedData.sigmax, fittedData.sigmay, fittedData.x0, fittedData.y0, X, Y);
fitDataTrunc = reshape(fitDataTrunc,size(subImg));
resid = fitDataTrunc - subImg;

figure
surf(resid)

%% Segmenting overlapping objects
clear
close all

% First we'll create two overlapping circles, centralized in opposing
% quadrants of an image
center1 = -10;
center2 = -center1;

% Make the radius large enough to overlap
dist = sqrt(2*(2*center1)^2);
radius = dist/2*1.4;

% Set the limits as the size of the image
lims = [floor(center1 - 1.2*radius) ceil(center2 + 1.2*radius)];

% Create a meshgrid to label our pixels
[x,y] = meshgrid(lims(1):lims(2));

% Define the two circles by sqrt(x^2 + y^2)
bw1 = sqrt((x-center1).^2 + (y - center1).^2) <= radius;
bw2 = sqrt((x-center2).^2 + (y - center2).^2) <= radius;

% Overlap the two circles
bw = bw1 | bw2;

figure
imshow(bw)

%%
% The circles are clearly overlapping, but represent two objects.  How can
% we tell Matlab that we want to ID this as two objects.  We can use the
% watershed algorithm.  Before we get there, we need to make the centers of
% the two objects the minima.  Then as the water fills in, it will start
% from there and fill in the objects.

% The bwdist image will calculate the distance to the neare
distImg = mat2gray(bwdist(~bw));

% Then we can take the complement image, and the centers of the circle are
% now the lowest points.  They will fill up with the watershed.
distImg2 = imcomplement(distImg);
figure
imshow([distImg distImg2])

% Now we need to set the background pixels back to zero, and we can run the
% watershed
distImg3 = distImg2.*bw;
figure
imshow(distImg3)
% Now we can call the watershed
L = watershed(distImg3);

figure
imshow(L,[])

%% Watershed on more real images
clear
close all

img = im2double(imread('QD_AFM.png'));

figure
imshow(img)

% Convert to gray image and threshold
img2 = rgb2gray(img);
img3 = mat2gray(imtophat(img2,strel('disk',2)));
img4 = imdilate(img3,strel('square',2));

bw = im2bw(img4,graythresh(img3));

figure
imshowpair(img2,bw)


distImg = mat2gray(bwdist(~bw));

% Then we can take the complement image, and the centers of the circle are
% now the lowest points.  They will fill up with the watershed.
distImg2 = imcomplement(distImg);
figure
imshow([distImg distImg2])

% Now we need to set the background pixels back to zero, and we can run the
% watershed
distImg3 = distImg2.*bw;
figure
imshow(distImg3)
% Now we can call the watershed
L = watershed(distImg3);
labelImg = label2rgb(L, @spring, 'c', 'shuffle');

figure
subplot(2,1,1)
imshow(labelImg,[])
subplot(2,1,2)
imshow(img)

%%
clear
close all

% Outline of how to fit STORM data
% 1. Create a high resolution single image for adding in spots, easy to
% make a 1 nm x 1 nm resolution

imgSTORM = zeros(ysize*100,xsize*100); % This is assuming 100 nm pixels

% 2. Load in the frame
for f = 1:nframes
    currImg = imread(file(:,:,j));
    

% 3. Reduce noise by convolving with a Gaussian with the known PSF width
    gFilt = fspecial('Gaussian',7,2);
    currImgFilt = conv2(imgTruncNoise,gFilt,'same');

    
% 4. Threshold the image to find your spots
    currImgBW = im2bw(currImgFilt,bwthresh(currImgFilt));
    
% 5. Use regionprops to label the spots
    stats = regionprops(currImgBW,currImg,'BoundingBox','MeanIntenisty','Area');
    bbs = [stats.BoundBox];
    avgIntens = [stats.MeanIntensity];
    areas = [stats.Area];
    nspots = length(stats);
    
% 6. Now lets cycle through the spots and fit each one to a Gaussian
    for s = 1:nspots
        subImg = currImg(boundbox_limits);
        
        %Initial guesses
        a1 = max(subImg(:));
        sigmax = 2;
        sigmay = 2;
        x0 = length(xPixels)/2;
        y0 = length(yPixels)/2;
        fitData = fit_2D_Gauss(subImg,a1,sigmax,sigmay,x0,y0);
        
 % 7. Calculate uncertainty width (new width for the Gaussian 
        fitUncer = (1/numPhotons)*(PSF^2 + pixSize^2/12) + ...
            (1/(pixSize*numPhotons^2))*(4*sqrt(pi)*PSF^3*b^2);
        
 % 8. Create gaussian spot 
        fitSpot = fspecial('Gaussian',100,fitUncer); % This will be a 100 nm x 100 nm Gaussian box
        
 % 9. Add this spot to the STORM image
        imgSTORM(yFit-50:yFit+49,xFit-50:xFit+49) = fitSpot;
    end
end

