intFindFociCurve

 intFindFociCurve: finds the location of foci and assigns them to the cells.
 The foci are found using the curve (L-) filtered image. The foci are the
 identified using watershed. The foci regions are assigned to the closest cell
 and the location of foci is found in subpixel resolution, using a gaussian fit.

 INPUT : 
       data : cell/regions file (err file)
       CONST : segmentation constants
       channelID : fluorescence channel number

 OUTPUT : 
       data : updated data with sub-pixel fluorescence model
       a locus[channelID] field is added with the following structure :
              r : sub-pixel global coordinates for foci location (x y)
              fitSigma : sigma of gaussian fit
              intensity : maximum intensity of gaussian fit
              fitScore : score of gaussian fitting to the foci (0 - 1)
              normIntensity : normalized intensity (divided by std(cell
              fluor)
              score : locus intensity / std(cell fluorescence) *  gaussian fit score;
              shortaxis : sub-pixel local coordinates for foci location (x y)
              longaxis : sub-pixel local coordinates for foci location (x y)

 Copyright (C) 2016 Wiggins Lab 
 Written by Connor Brennan, Stella Stylianidou & Paul Wiggins.
 University of Washington, 2016
 This file is part of SuperSegger.
 
 SuperSegger is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 SuperSegger is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with SuperSegger.  If not, see <http://www.gnu.org/licenses/>.

0001 function [ data ] = intFindFociCurve( data, CONST, channelID )
0002 % intFindFociCurve: finds the location of foci and assigns them to the cells.
0003 % The foci are found using the curve (L-) filtered image. The foci are the
0004 % identified using watershed. The foci regions are assigned to the closest cell
0005 % and the location of foci is found in subpixel resolution, using a gaussian fit.
0006 %
0007 % INPUT :
0008 %       data : cell/regions file (err file)
0009 %       CONST : segmentation constants
0010 %       channelID : fluorescence channel number
0011 %
0012 % OUTPUT :
0013 %       data : updated data with sub-pixel fluorescence model
0014 %       a locus[channelID] field is added with the following structure :
0015 %              r : sub-pixel global coordinates for foci location (x y)
0016 %              fitSigma : sigma of gaussian fit
0017 %              intensity : maximum intensity of gaussian fit
0018 %              fitScore : score of gaussian fitting to the foci (0 - 1)
0019 %              normIntensity : normalized intensity (divided by std(cell
0020 %              fluor)
0021 %              score : locus intensity / std(cell fluorescence) *  gaussian fit score;
0022 %              shortaxis : sub-pixel local coordinates for foci location (x y)
0023 %              longaxis : sub-pixel local coordinates for foci location (x y)
0024 %
0025 % Copyright (C) 2016 Wiggins Lab
0026 % Written by Connor Brennan, Stella Stylianidou & Paul Wiggins.
0027 % University of Washington, 2016
0028 % This file is part of SuperSegger.
0029 %
0030 % SuperSegger is free software: you can redistribute it and/or modify
0031 % it under the terms of the GNU General Public License as published by
0032 % the Free Software Foundation, either version 3 of the License, or
0033 % (at your option) any later version.
0034 %
0035 % SuperSegger is distributed in the hope that it will be useful,
0036 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0037 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0038 % GNU General Public License for more details.
0039 %
0040 % You should have received a copy of the GNU General Public License
0041 % along with SuperSegger.  If not, see <http://www.gnu.org/licenses/>.
0042 
0043 DEBUG_FLAG = false;
0044 fieldname = ['locus', num2str(channelID)];
0045 options =  optimset('MaxIter', 1000, 'Display', 'off', 'TolX', 1/10);
0046 
0047 % Get images out of the structures.
0048 originalImage = double(data.(['fluor',num2str(channelID)]));
0049 
0050 % Subtract gaussian blurred image
0051 % to get rid of big structure
0052 hg = fspecial( 'gaussian' , 210, 30 );
0053 highPassImage = originalImage - imfilter( originalImage, hg, 'replicate' );
0054 cytoplasmicFlourescenceSTD = std(double(highPassImage(data.mask_bg)));
0055 if isnan(cytoplasmicFlourescenceSTD)
0056     cytoplasmicFlourescenceSTD = 1;
0057     disp ('Possibly empty mask for image - foci may not be found');
0058 end
0059 normalizedImage = originalImage/cytoplasmicFlourescenceSTD; % normalization so that intensities;
0060 
0061 %Take only pixels above 1 std (noise reduction?)
0062 normalizedImage = normalizedImage - 1;
0063 normalizedImage(normalizedImage < 0) = 0; % logical mask of foci found
0064 
0065 filter_width = 1.5;
0066 [~,~,fluorFiltered] = curveFilter(normalizedImage,filter_width);
0067 data.(['fluor',num2str(channelID),'_filtered']) = fluorFiltered;
0068 
0069 mask_mod = bwmorph (data.mask_bg, 'dilate', 1);
0070 
0071 fociWatershed = watershed(-fluorFiltered); % watershed to identify foci
0072 maskedFociWatershed = logical(double(fociWatershed).*double(mask_mod));
0073 
0074 fociRegionLabels = bwlabel(maskedFociWatershed);
0075 props = regionprops( fociRegionLabels, {'BoundingBox'} );
0076 numFociRegions = numel(props);
0077 imsize = size(highPassImage);
0078 
0079 % initialize focus fields
0080 focusInit.r = [nan,nan];
0081 focusInit.score = nan;
0082 focusInit.intensity = nan;
0083 focusInit.normIntensity = nan;
0084 focusInit.shortaxis = nan;
0085 focusInit.longaxis = nan;
0086 focusInit.fitSigma = nan;
0087 focusInit.fitScore = nan;
0088 
0089 fociData = [];
0090 cellIDs = [];
0091 
0092 if DEBUG_FLAG
0093     figure(2);
0094     clf;
0095 end
0096 
0097 for ii = 1:numFociRegions
0098     tempData = focusInit;
0099     
0100     [xPad,yPad] = getBBpad( props(ii).BoundingBox, imsize, 3 );
0101     [meshX,meshY] = meshgrid(xPad, yPad);
0102     
0103     maskToFit = (fociRegionLabels(yPad, xPad) == ii); % foci region
0104     imageToFit  = fluorFiltered(yPad, xPad);  % filtered image
0105     imageToFit = imageToFit .* double(maskToFit);
0106     
0107     [~, maxIndex] = max(imageToFit(maskToFit));
0108     
0109     tempImage = imageToFit(maskToFit);
0110     tempData.intensity = tempImage(maxIndex);
0111    
0112     tempImage = meshX(maskToFit);
0113     fociX = tempImage(maxIndex);
0114     tempData.r(1) = fociX;
0115    
0116     tempImage = meshY(maskToFit);
0117     fociY = tempImage(maxIndex);
0118     tempData.r(2) = fociY;
0119 
0120     % figure out which cell the focus belongs to
0121     maskSize = [numel(yPad),numel(xPad)];
0122 
0123    
0124     cellsLabel = data.regs.regs_label(yPad,xPad);
0125     cellsMask = logical(cellsLabel);
0126     tempMask = zeros(maskSize);
0127     
0128     tempMask(fociY - yPad(1)+1, fociX - xPad(1)+1 ) = 1;
0129     distanceToFoci = bwdist( tempMask );
0130 
0131     cellIDList = cellsLabel(cellsMask);
0132     [~, minDistanceIndex] = min(distanceToFoci(cellsMask));
0133     bestCellID = cellIDList(minDistanceIndex);
0134     
0135     if ~isempty( bestCellID )
0136         croppedImage = highPassImage(data.CellA{bestCellID}.yy, data.CellA{bestCellID}.xx);
0137         cellFlouresenseSTD = std(croppedImage(data.CellA{bestCellID}.mask));  
0138         
0139         if tempData.intensity >0
0140             %Initialize parameters
0141             backgroundIntensity = 0;
0142             gaussianIntensity = fluorFiltered(fociY, fociX) - backgroundIntensity;
0143             sigmaValue = 1;
0144 
0145             parameters(1) = fociX;
0146             parameters(2) = fociY;
0147             parameters(3) = gaussianIntensity;
0148             parameters(4) = sigmaValue;
0149             %parameters(5) = backgroundIntensity;
0150 
0151             [parameters] = fminsearch( @doFit, parameters, options);
0152 
0153             gaussianApproximation = makeGassianTestImage(meshX, meshY, parameters(1), parameters(2), parameters(3), backgroundIntensity, parameters(4));
0154 
0155             %Crop out fit gaussian from original image
0156             croppedImage = imageToFit;
0157             croppedImage(gaussianApproximation < 0.1 * max(max(gaussianApproximation))) = 0;
0158 
0159             imageTotal = sqrt(sum(sum(croppedImage)));
0160             guassianTotal = sqrt(sum(sum(gaussianApproximation)));
0161 
0162             fitScore = sum(sum(sqrt(croppedImage) .* sqrt(gaussianApproximation))) / (imageTotal * guassianTotal);
0163 
0164             if DEBUG_FLAG
0165                 figure(1);
0166                 clf;
0167                 subplot(2, 2, 1);
0168                 imshow(imageToFit, []);
0169                 subplot(2, 2, 2);        
0170                 imshow(croppedImage, []);    
0171                 subplot(2, 2, 3);        
0172                 imshow(gaussianApproximation, []);   
0173                 subplot(2, 2, 4);          
0174                 title(['Score: ', num2str(fitScore)]);
0175                 keyboard;
0176             end
0177 
0178             tempData.r(1) = parameters(1);
0179             tempData.r(2) = parameters(2);
0180             tempData.fitSigma = parameters(4);
0181             tempData.intensity = parameters(3);
0182             tempData.fitScore = fitScore;
0183 
0184             %Calculate scores
0185             tempData.normIntensity = tempData.intensity;
0186             tempData.score = tempData.intensity / (sqrt(pi)*tempData.fitSigma) * tempData.fitScore;
0187             tempData.shortaxis = ...
0188                 (tempData.r-data.CellA{bestCellID}.coord.rcm)*data.CellA{bestCellID}.coord.e2;
0189             tempData.longaxis = ...
0190                 (tempData.r-data.CellA{bestCellID}.coord.rcm)*data.CellA{bestCellID}.coord.e1;
0191 
0192 
0193             %Assign to array
0194             cellIDs(ii) = bestCellID;  
0195             focusData(ii) = tempData;
0196             if DEBUG_FLAG
0197                figure(2);
0198                hold on;
0199                plot(fociX, fociY, '.r' );
0200                text(fociX, fociY, num2str( tempData.intensity, '%1.2g' ));
0201             end
0202         end
0203     end
0204 end
0205 
0206 
0207 % assign to cells
0208 for ii = 1:data.regs.num_regs
0209 
0210     fociIndex = find(cellIDs == ii);
0211     tempData = focusData(fociIndex);
0212     
0213     [~, order] = sort( [tempData.intensity], 'descend' );
0214     sortedFoci = tempData(order);
0215     
0216     focus = focusInit;
0217     if numel(sortedFoci) > 0
0218         maxIndex = find([sortedFoci.intensity] > 0.333 * sortedFoci(1).intensity);
0219         if numel(maxIndex) > CONST.trackLoci.numSpots(channelID)
0220             maxIndex = maxIndex(1:CONST.trackLoci.numSpots(channelID));
0221         end
0222 
0223         numFoci = numel(maxIndex);
0224 
0225         for jj = 1:numFoci
0226             focus(jj) = sortedFoci(jj);
0227         end 
0228     end
0229     
0230     scores = [focus(:).score];
0231     focus = focus(~isnan(scores));  
0232     data.CellA{ii}.(fieldname) = focus;
0233     
0234     % creates a filtered image of the cell
0235     xPad = data.CellA{ii}.xx;
0236     yPad = data.CellA{ii}.yy;
0237     data.CellA{ii}.(['fluor',num2str(channelID),'_filtered'])=fluorFiltered( yPad, xPad );
0238 end
0239 
0240     function error = doFit(parameters)
0241         % doFit : does the gaussian fit to the foci and calculates the
0242         % error.
0243         %parameters are :
0244         %parameter(1) - Sub-pixel resolution of foci position X
0245         %parameter(2) - Sub-pixel resolution of foci position Y
0246         %parameter(3) - Intensity of the gaussian
0247         %parameter(4) - sigma of gaussian
0248         %parameter(5) - background intensity
0249         gaussian = makeGassianTestImage(meshX, meshY, parameters(1), parameters(2), parameters(3), backgroundIntensity, parameters(4));
0250         tempImage = (double(imageToFit) - gaussian);
0251         error = sum(sum(tempImage.^2));
0252     end
0253 
0254     function testImage = makeGassianTestImage(meshX, meshY, fociX, fociY, gaussianIntensity, backgroundIntensity, sigmaValue)
0255         testImage = backgroundIntensity + gaussianIntensity * exp( -((meshX - fociX).^2 + (meshY - fociY).^2)/(2 * sigmaValue^2) );
0256     end
0257 
0258 end