Description of cellprops3

0001 function info = cellprops3( mask, props )
0002 % cellprops3 : Calculates the shape properties of the region or 'cell'.
0003 %
0004 % INPUT :
0005 %    mask : cell / region mask
0006 %    props : properties generated by regionprops
0007 % OUTPUT :
0008 %   info = { 'long axis: L1', ...
0009 %     'short axis mean:  L2mean', ...
0010 %     'neck width:  Lneck,
0011 %     'short axis max: L2max,
0012 %     'short axis var: L2v', ...
0013 %     '<width sin pi s/S >: s1 ', ...
0014 %     '<width sin pi 2s/S>: s2', ...
0015 %     '<width sin pi 3s/S>: s3', ...
0016 %     'RoundIndOver: RoundIndOver', ...
0017 %     'RoundIndUnder RoundIndUnder', ...
0018 %     'Area: props.Area,', ....
0019 %     'max dtheta^2: stm1 ', ...
0020 %     '<theta sin pi s/2S >: sa1 ', ...
0021 %     '<theta sin pi s/S  >: sa2', ...
0022 %     '<theta sin pi 3s/2S>: sa3', ...
0023 %     '<theta sin pi 2s/S >: sa4', ...
0024 %     '1/long axis: 1/L1', ...
0025 %     '1/mean short axis: 1/L2mean',...
0026 %     'Min End Width: minWidthEnd',...
0027 %     'round end min: sqmin',...
0028 %     'round end max: sqmax'};
0029 %
0030 % Copyright (C) 2016 Wiggins Lab
0031 % Written by Stella Stylianidou, Paul Wiggins.
0032 % University of Washington, 2016
0033 % This file is part of SuperSegger.
0034 %
0035 % SuperSegger is free software: you can redistribute it and/or modify
0036 % it under the terms of the GNU General Public License as published by
0037 % the Free Software Foundation, either version 3 of the License, or
0038 % (at your option) any later version.
0039 %
0040 % SuperSegger is distributed in the hope that it will be useful,
0041 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0042 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0043 % GNU General Public License for more details.
0044 %
0045 % You should have received a copy of the GNU General Public License
0046 % along with SuperSegger.  If not, see <http://www.gnu.org/licenses/>.
0047 
0048 debug = false;
0049 
0050 
0051 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0052 %
0053 % Calculate cell thickness and length by rotating the
0054 % region based on the angle of the majorAxis.
0055 %
0056 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0057 Orientation = props.Orientation;
0058 
0059 % rotated image
0060 imRot = fast_rotate_loose_double( mask, -Orientation+90);
0061 ss = size(imRot);
0062 
0063 if debug
0064     figure(2);
0065     clf;
0066     imshow(cat(3,ag(imRot)+0.5*ag(~imRot),ag(imRot),ag(imRot)),'InitialMagnification','fit');
0067     hold on;
0068 end
0069 
0070 % short axis width
0071 width       = sum(double(imRot),2);
0072 ymask       = (width>=1);
0073 widthWindow = ([0;width(1:end-1)] + width + [width(2:end);0] )/3;
0074 
0075 L1          = sum(ymask); % long axis
0076 L2max       = max(widthWindow); % max of short axis
0077 L2v         = var(widthWindow); % variance of short axis
0078 
0079 ymin = find(ymask,1,'first');
0080 ymax = find(ymask,1,'last');
0081 ind = (ymin):(ymax);
0082 L2mean = mean(width(ind)); % mean of short axis
0083 
0084 if L2mean<1
0085     L2mean=1;
0086 end
0087 
0088 
0089 % sine cosine idea : width(0:L1) * sin (2*pi*(0:L1)/L1)
0090 
0091 nn = (numel(ind)-1); % long axis width
0092 s  = 2*pi*(0:nn)/nn;
0093 s1 = mean(width(ind)'.*sin(s/2));
0094 s2 = mean(width(ind)'.*sin(s));
0095 s3 = mean(width(ind)'.*sin(s/3));
0096 
0097 
0098 % neck width
0099 if numel(widthWindow) > 6;
0100     ww       = widthWindow';
0101     v        = ww(3:end)-ww(1:end-2); % change in width from i3 to i1
0102     
0103     s_change = logical([0, sign(v(2:end))-sign(v(1:end-1)),0,0]); % positive change in change
0104     s_change(2:end) = logical(s_change(2:end)+s_change(1:end-1));
0105     ind = find(s_change);
0106     ind_ = ind(logical((ind>3).*(ind<(numel(ww)-2))));
0107     
0108     if isempty(ind_)
0109         Lneck = L2mean;
0110         pos = mean([ymin,ymax]);
0111     else
0112         [Lneck,pos] = min(ww(ind_)); % minimum width value
0113         pos = ind_(pos); % index in width window
0114     end
0115 else
0116         Lneck = L2mean;
0117         pos = mean([ymin,ymax]);
0118 end
0119 
0120 if Lneck > L2mean
0121     Lneck = L2mean;
0122     pos = mean([ymin,ymax]);
0123 end
0124 
0125 if debug
0126     % ww : width window
0127     xx = 1:numel(ww);
0128     figure(1);
0129     clf;
0130     plot( xx, ww, '.-' );
0131     ylabel ('width');
0132     xlabel ('Image Pixel (1 : width window)');
0133     hold on;
0134     x = 1:Lneck
0135     y=pos*ones(1,numel(x))
0136     plot(y,x, 'r--'); % marks position of neck
0137     plot( xx(ind_), ww(ind_), 'g.' );
0138     figure(2)
0139     
0140 end
0141 
0142 
0143 % finds thin ends
0144 im = find( logical(width), 1, 'first' );
0145 ip = find( logical(width), 1, 'last' );
0146 if ip-im >= 3
0147    m1 = mean( width(im+[ 0, 1]) );
0148    m4 = mean( width(ip+[-1, 0]) );
0149    minWidthEnd = min(m1,m4);
0150 else
0151    minWidthEnd = L2mean;
0152 end
0153 
0154 % finds square ends
0155 if ip-im >= 7
0156     m2 = mean( width(im+[ 2, 3]) );
0157     m3 = mean( width(ip+[-3,-2]) );
0158 
0159     sq1 = m2-m1;
0160     sq2 = m3-m4;
0161     
0162     sqmax = max(sq1,sq2);
0163     sqmin = min(sq1,sq2);
0164 else
0165    sqmax = 0;
0166    sqmin = 0;
0167 end
0168 
0169 % maximum bend angle
0170 dr  = 3;
0171 dr2 = 3;
0172 yy = (ymin+dr):dr2:(ymax-dr);
0173 
0174 if numel(yy) > 4
0175     tmp = imRot(yy,:);
0176     
0177     x = 1:ss(2);
0178     
0179     [X,Y] = meshgrid(x,yy);
0180     
0181     xx = sum( X.*tmp,2 )./sum(tmp,2);
0182     
0183     d1 = [xx(1:end-2)-xx(3:end),yy(1:end-2)'-yy(3:end)'];
0184     
0185     dangle = ((d1(1:end-2,1).*d1(3:end,2)...
0186         -d1(1:end-2,2).*d1(3:end,1))./...
0187         (norm2(d1(1:end-2,:)).*norm2(d1(3:end,:))));
0188             
0189     stm1 = max( dangle.^2 );
0190     
0191     
0192     angle = cumsum(dangle)';
0193     na = numel(angle)-1;
0194     
0195     s = (0:na)/na*pi/2;
0196 
0197     try
0198     sa1 = mean( angle.*sin(s  ) );
0199     sa2 = mean( angle.*sin(s*2) );
0200     sa3 = mean( angle.*sin(s*3) );
0201     sa4 = mean( angle.*sin(s*4) );
0202     catch ME
0203         printError(ME);
0204     end     
0205         
0206     if debug       
0207         ss_ = size(d1);       
0208         for ii = 1:ss_(1)            
0209             plot( xx(ii)+[0,-d1(ii,1)], yy(ii)+[0,-d1(ii,2)], 'c.' );
0210         end
0211     end
0212     
0213     
0214     
0215     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0216     %
0217     % Calculate the end score
0218     %
0219     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0220     yy = [(ymin+dr),(ymax-dr)];
0221     xx = widthWindow(round(yy));
0222     r = ((xx)/2);
0223     r(r<1)=1;
0224     
0225     yyd = [(ymin+r(1)+[-1,0,1]),(ymax-r(2)+[-1,0,1])];
0226     yy  = round(yyd);
0227     
0228     [X] = meshgrid(x,yy);
0229     
0230     tmp = imRot(yy,:);
0231     xx = sum( X.*tmp,2 )./sum(tmp,2);
0232     xx = [mean(xx(1:3)),mean(xx(4:6))];
0233     yyd = yyd([2,5]);
0234 
0235     
0236     if debug
0237         phi = 0:.1:2*pi;
0238         
0239         plot(xx,yyd,'go');
0240         xxx = r(1)*cos(phi);
0241         yyy = r(1)*sin(phi);
0242         plot(xx(1)+xxx,yy(1)+yyy,'g:')
0243         
0244         xxx = r(2)*cos(phi);
0245         yyy = r(2)*sin(phi);
0246         plot(xx(2)+xxx,yyd(2)+yyy,'g:')
0247         
0248         plot( [1,ss(2)],[ymax,ymax],':r')
0249         plot( [1,ss(2)],[ymin,ymin],':r')
0250         
0251     end
0252     
0253     xcen = xx;
0254     ycen = yyd;
0255     
0256     x0 = xcen(1);
0257     y0 = ycen(1);
0258     r0 = r(1);
0259     
0260     yy = max(1,round(y0-r0-2)):min(round(y0+r0+2),ss(1));
0261     xx = max(1,round(x0-r0-2)):min(round(x0+r0+2),ss(2));
0262     
0263     [X,Y] = meshgrid(xx,yy);
0264     tmp = imRot(yy,xx);
0265     
0266     R2 = (X-x0).^2 + (Y-y0).^2 - r0^2;
0267     
0268     ind_over = find( ((R2 > 0).*tmp).*(d1(1,1)*(X-x0) + d1(1,2)*(Y-y0) > 0));
0269     ind_under = find( ((R2 < 0).*(tmp<1)).*(d1(1,1)*(X-x0) + d1(1,2)*(Y-y0) > 0));
0270     
0271     
0272     if debug
0273         plot(X(ind_over),Y(ind_over),'r.');
0274         plot(X(ind_under),Y(ind_under),'b.');
0275     end
0276     
0277     RoundIndOverTop  = sum( tmp(ind_over).*R2(ind_over))/r0^4;
0278     RoundIndUnderTop = sum( (1-tmp(ind_under)).*abs(R2(ind_under)))/r0^4;
0279     
0280     x0 = xcen(2);
0281     y0 = ycen(2);
0282     r0 = r(2);
0283     
0284     yy = max(1,round(y0-r0-2)):min(round(y0+r0+2),ss(1));
0285     xx = max(1,round(x0-r0-2)):min(round(x0+r0+2),ss(2));
0286     
0287     [X,Y] = meshgrid(xx,yy);
0288     tmp = imRot(yy,xx);
0289     
0290     R2 = (X-x0).^2 + (Y-y0).^2 - r0^2;
0291     
0292     ind_over = find( ((R2 > 0).*tmp).*(d1(1,1)*(X-x0) + d1(1,2)*(Y-y0) < 0));
0293     ind_under = find( ((R2 < 0).*(tmp<1)).*(d1(1,1)*(X-x0) + d1(1,2)*(Y-y0) < 0));
0294     
0295     if debug
0296         plot(X(ind_over),Y(ind_over),'r.');
0297         plot(X(ind_under),Y(ind_under),'b.');
0298     end
0299     
0300     
0301     RoundIndOverBot  = sum( tmp(ind_over).*R2(ind_over))/r0^4;
0302     RoundIndUnderBot = sum( (1-tmp(ind_under)).*abs(R2(ind_under)))/r0^4;
0303     
0304 else
0305     %'fail'
0306     stm1 = 0;
0307     sa1 = 0;
0308     sa2 = 0;
0309     sa3 = 0;
0310     sa4 = 0;
0311     RoundIndOverTop  = 0;
0312     RoundIndUnderTop = 0;
0313     RoundIndOverBot  = 0;
0314     RoundIndUnderBot = 0;
0315     
0316 end
0317 
0318 
0319 RoundIndOver  = RoundIndOverTop  + RoundIndOverBot;
0320 RoundIndUnder = RoundIndUnderTop + RoundIndUnderBot;
0321 
0322 
0323 
0324 
0325 if L1 < 1
0326     L1 = 1;
0327 end
0328 
0329 
0330 info = [L1, ...
0331     L2mean, ...
0332     Lneck, ...
0333     L2max, ...
0334     L2v, ...
0335     s1, ...
0336     s2, ...
0337     s3, ...
0338     RoundIndOver, ...
0339     RoundIndUnder, ...
0340     props.Area, ....
0341     stm1, ...
0342     sa1, ...
0343     sa2, ...
0344     sa3, ...
0345     sa4, ...
0346     1/L1, ...
0347     1/L2mean,...
0348     minWidthEnd,...
0349     sqmin,...
0350     sqmax];
0351 
0352 info(isnan(info)) = 0;
0353 
0354 if numel(info) ~= 21
0355     disp ('incorrect number of cell properties calculated')
0356 end
0357 
0358 
0359 
0360 if debug
0361     drawnow;
0362     keyboard;
0363 end
0364 
0365 
0366 
0367 end
0368 
0369 function tmp = norm2( dd )
0370 
0371 tmp = sqrt(sum(dd.*dd,2));
0372 
0373 end
cellprops3

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
