Description of CellsortSegmentation

0001 function [ica_segments, segmentlabel, segcentroid] = CellsortSegmentation(ica_filters, smwidth, thresh, arealims, plotting)
0002 % [ica_segments, segmentlabel, segcentroid] = CellsortSegmentation(ica_filters, smwidth, thresh, arealims, plotting)
0003 %
0004 %CellsortSegmentation
0005 % Segment spatial filters derived by ICA
0006 %
0007 % Inputs:
0008 %     ica_filters - X x Y x nIC matrix of ICA spatial filters
0009 %     smwidth - standard deviation of Gaussian smoothing kernel (pixels)
0010 %     thresh - threshold for spatial filters (standard deviations)
0011 %     arealims - 2-element vector specifying the minimum and maximum area
0012 %     (in pixels) of segments to be retained; if only one element is
0013 %     specified, use this as the minimum area
0014 %     plotting - [0,1] whether or not to show filters
0015 %
0016 % Outputs:
0017 %     ica_segments - segmented spatial filters
0018 %     segmentabel - indices of the ICA filters from which each segment was derived
0019 %     segcentroid - X,Y centroid, in pixels, of each segment
0020 %
0021 % Eran Mukamel, Axel Nimmerjahn and Mark Schnitzer, 2009
0022 % Email: eran@post.harvard.edu, mschnitz@stanford.edu
0023 
0024 tic
0025 fprintf('-------------- CellsortSegmentation %s -------------- \n', date)
0026 
0027 if (nargin<3)||isempty(thresh)
0028     thresh = 2;
0029 end
0030 if (nargin<4)||isempty(arealims)
0031     arealims = 200;
0032 end
0033 if (nargin<5)||isempty(plotting)
0034     plotting = 0;
0035 end
0036 [nic,pixw,pixh] = size(ica_filters);
0037 
0038 ica_filtersorig = ica_filters / abs(std(ica_filters(:)));
0039 ica_filters = (ica_filters - mean(ica_filters(:)))/abs(std(ica_filters(:)));
0040 if smwidth>0
0041     % Smooth mixing filter with a Gaussian of s.d. smwidth pixels
0042     smrange = max(5,3*smwidth);
0043     [x,y] = meshgrid([-smrange:smrange]);
0044 
0045     smy = 1; smx = 1;
0046     ica_filtersfilt = exp(-((x/smx).^2 + (y/smy).^2)/(2*smwidth^2));
0047     
0048     ica_filtersfilt = ica_filtersfilt/sum(ica_filtersfilt(:));
0049     ica_filtersbw = false(pixw,pixh,nic);
0050     tic
0051     for j = 1:size(ica_filters,1)
0052         ica_filtersuse = ica_filters(j,:,:);
0053         ica_filtersuse = (ica_filtersuse - mean(ica_filtersuse(:)))/abs(std(ica_filtersuse(:)));
0054         ica_filtersbw(:,:,j) = (imfilter(ica_filtersuse, ica_filtersfilt, 'replicate', 'same') > thresh);
0055     end
0056 else
0057     ica_filtersbw = (permute(ica_filters,[2,3,1]) > thresh);
0058     ica_filtersfilt = 1;
0059 end
0060 
0061 tic
0062 if plotting
0063     clf
0064     set(gcf,'Color','w')
0065     colormap(gray)
0066     
0067     subplot(223)
0068     imagesc(squeeze(sum(ica_filters,1)))
0069     axis image off
0070     hold on
0071 end
0072 ica_filterslabel = [];
0073 ica_segments = [];
0074 k=0;
0075 L=[];
0076 segmentlabel = [];
0077 segcentroid = [];
0078 [x,y] = meshgrid([1:pixh], [1:pixw]);
0079 for j = 1:nic
0080     % Label contiguous components
0081     L = bwlabel(ica_filtersbw(:,:,j), 4);
0082     Lu = 1:max(L(:));
0083 
0084     % Delete small components
0085     Larea = struct2array(regionprops(L, 'area'));
0086     Lcent = regionprops(L, 'Centroid');
0087     
0088     if length(arealims)==2
0089         Lbig = Lu( (Larea >= arealims(1))&(Larea <= arealims(2)));
0090         Lsmall = Lu((Larea < arealims(1))|(Larea > arealims(2)));
0091     else
0092         Lbig = Lu(Larea >= arealims(1));
0093         Lsmall = Lu(Larea < arealims(1));
0094     end
0095     
0096     L(ismember(L,Lsmall)) = 0;   
0097         
0098     for jj = 1:length(Lbig)
0099         segcentroid(jj+k,:) = Lcent(Lbig(jj)).Centroid;
0100     end
0101     
0102     ica_filtersuse = squeeze(ica_filtersorig(j,:,:));
0103     for jj = 1:length(Lbig)
0104         ica_segments(jj+k,:,:) = ica_filtersuse .* ( 0*(L==0) + (L==Lbig(jj)) );  % Exclude background
0105     end
0106     
0107     if plotting && ~isempty(Lbig)
0108         if smwidth>0
0109             subplot(2,2,2)
0110             ica_filtersuse = squeeze(ica_filters(j,:,:));
0111             ica_filtersuse = (ica_filtersuse - mean(ica_filtersuse(:)))/abs(std(ica_filtersuse(:)));
0112             imagesc(imfilter((ica_filtersuse), ica_filtersfilt, 'replicate', 'same'),[-1,4])
0113             hold on
0114             contour(imfilter((ica_filtersuse), ica_filtersfilt, 'replicate', 'same'), [1,1]*thresh, 'k')
0115             hold off
0116             hc = colorbar('Position',[0.9189    0.6331    0.0331    0.2253]);
0117             ylabel(hc,'Std. dev.')
0118             title(['IC ',num2str(j),' smoothed'])
0119             axis image off
0120 
0121             subplot(2,2,1)
0122         else
0123             subplot(211)
0124         end
0125         imagesc(squeeze(ica_filters(j,:,:)))
0126         title(['IC ',num2str(j),' original'])
0127         axis image off
0128 
0129         colord = lines(k+length(Lbig));
0130         for jj = 1:length(Lbig)
0131             subplot(223)
0132             contour(ica_filtersbw(:,:,j), [1,1]*0.5, 'color',colord(jj+k,:),'linewidth',2)
0133             hold on
0134             text(segcentroid(jj+k,1), segcentroid(jj+k,2), num2str(jj+k), 'horizontalalignment','c', 'verticalalignment','m')
0135             set(gca, 'ydir','reverse','tickdir','out')
0136             axis image
0137             xlim([0,pixw]); ylim([0,pixh])
0138             
0139             subplot(224)
0140             imagesc(squeeze(ica_segments(jj+k,:,:)))
0141             hold on
0142             plot(segcentroid(jj+k,1), segcentroid(jj+k,2), 'bo')
0143             hold off
0144             axis image off
0145             title(['Segment ',num2str(jj+k)])
0146             drawnow
0147         end
0148     end
0149     k = size(ica_segments,1);
0150 end
0151 toc
CellsortSegmentation

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
