cellprops5

 cellprops5 : Calculates the shape properties of a region or 'cell'.
 This are used later to calculate the score of a region.

 INPUT :
       mask : cell / region mask
       props : properties generated by regionprops

0001 function info = cellprops5( mask, props )
0002 % cellprops5 : Calculates the shape properties of a region or 'cell'.
0003 % This are used later to calculate the score of a region.
0004 %
0005 % INPUT :
0006 %       mask : cell / region mask
0007 %       props : properties generated by regionprops
0008 
0009 % OUTPUT :
0010 %         info = [L1 : long axis of the region
0011 %             L2mean : mean of the short axis
0012 %             Lneck : neck width
0013 %             L2max : short axis max
0014 %             L2v : variance of the short axis
0015 %             s1 : <width sin pi s/S >
0016 %             s2 : <width sin pi 2s/S >
0017 %             s3 : <width sin pi 3s/S >
0018 %             RoundIndOver :
0019 %             RoundIndUnder,
0020 %             Area,
0021 %             stm1 : max dtheta^2
0022 %             sa1 : <theta sin pi s/2S >
0023 %             sa2 : <theta sin pi s/S >
0024 %             sa3 : <theta sin pi 3s/2S >
0025 %             sa4 : <theta sin pi 2s/S >
0026 %             1/L1 : 1 / long axis
0027 %             1/L2mean : 1 / mean of short axis
0028 %             minWidthEnd : Min End Width
0029 %             sqmin : round end min
0030 %             sqmax : round end max];
0031 %
0032 % Copyright (C) 2016 Wiggins Lab
0033 % University of Washington, 2016
0034 % This file is part of SuperSeggerOpti.
0035 
0036 debug = false;
0037 
0038 % Calculate cell thickness and length by rotating the
0039 % region based on the angle of the majorAxis.
0040 Orientation = props.Orientation;
0041 
0042 imRot = fast_rotate_loose_double( mask, -Orientation+90);
0043 ss = size(imRot);
0044 
0045 if debug
0046     figure(2);
0047     clf;
0048     imshow(cat(3,ag(imRot)+0.5*ag(~imRot),ag(imRot),ag(imRot)),'InitialMagnification','fit');
0049     hold on;
0050 end
0051 
0052 width = sum(double(imRot),2); % array of long axis length, per y direction
0053 ymask = (width>=1); % 0 if no cell along x-dimension, 1 if there is
0054 widthWindow = ([0;width(1:end-1)] + width + [width(2:end);0] )/3;
0055 % width in y-before + width in current y + width in y after /3
0056 
0057 L1          = sum(ymask);
0058 L2max       = max(widthWindow);
0059 L2v         = var(widthWindow);
0060 
0061 ymin = find(ymask,1,'first');
0062 ymax = find(ymask,1,'last');
0063 
0064 
0065 ind = (ymin):(ymax);
0066 L2mean = mean(width(ind));
0067 if L2mean<1
0068     L2mean=1;
0069 end
0070 
0071 % Sine Cosine Idea
0072 nn = (numel(ind)-1);
0073 s  = 2*pi*(0:nn)/nn;
0074 
0075 s1 = mean(width(ind)'.*sin(s*0.5));
0076 s2 = mean(width(ind)'.*sin(s));
0077 s3 = mean(width(ind)'.*sin(s*1.5));
0078 
0079 
0080 % Calculate the neck width
0081 if numel(widthWindow) > 6;
0082     ww       = widthWindow';
0083     v        = ww(3:end)-ww(1:end-2);
0084     s_change = logical([0, sign(v(2:end))-sign(v(1:end-1)),0,0]);
0085     s_change(2:end) = logical(s_change(2:end)+s_change(1:end-1));
0086     ind = find(s_change);
0087     ind_ = ind(logical((ind>3).*(ind<(numel(ww)-2))));
0088     
0089     if isempty(ind_)
0090         Lneck = L2mean;
0091         pos = mean([ymin,ymax]);
0092     else
0093         [Lneck,pos] = min( ww(ind_) );
0094         pos = ind_(pos);
0095     end
0096 else
0097     Lneck = L2mean;
0098     pos = mean([ymin,ymax]);
0099 end
0100 
0101 if Lneck > L2mean
0102     Lneck = L2mean;
0103     pos = mean([ymin,ymax]);
0104 end
0105 
0106 if debug
0107     xx = 1:numel(ww);    
0108     figure(1);
0109     clf;
0110     plot( xx, ww, '.-' );
0111     hold on;
0112     plot( pos+0*x,x, 'r:' );
0113     plot( xx(ind_), ww(ind_), 'g.' );
0114     figure(2)
0115 end
0116 
0117 % make stuff w/o rotating
0118 if abs(tan( props.Orientation*pi/180 )) > 1
0119     yunr = sum(mask,2)';
0120 else
0121     yunr = sum(mask);
0122 end
0123 
0124 ip = find( logical(yunr), 1, 'last' );
0125 im = find( logical(yunr), 1, 'first' );
0126 
0127 
0128 % thin ends
0129 im = find( logical(yunr), 1, 'first' );
0130 ip = find( logical(yunr), 1, 'last' );
0131 
0132 if ip-im >= 3
0133     m1 = mean( yunr(im+[ 0, 1]) );
0134     m4 = mean( yunr(ip+[-1, 0]) );
0135     minWidthEnd = min(m1,m4);
0136 else
0137     minWidthEnd = L2mean;
0138 end
0139 
0140 
0141 
0142 % Square ends
0143 if ip-im >= 7
0144     m2 = mean( yunr(im+[ 2, 3]) );
0145     m3 = mean( yunr(ip+[-3,-2]) );
0146     
0147     sq1 = m2-m1;
0148     sq2 = m3-m4;
0149     
0150     sqmax = max(sq1,sq2);
0151     sqmin = min(sq1,sq2);
0152 else
0153     sqmax = 0;
0154     sqmin = 0;
0155 end
0156 
0157 
0158 
0159 % Calculate the maximum bend angle.
0160 dr  = 3;
0161 dr2 = 3;
0162 yy = (ymin+dr):dr2:(ymax-dr);
0163 
0164 if numel(yy) > 4
0165     tmp = imRot(yy,:);
0166     x = 1:ss(2);    
0167     [X,Y] = meshgrid(x,yy);
0168     
0169     xx = sum( X.*tmp,2 )./sum(tmp,2);   
0170     d1 = [xx(1:end-2)-xx(3:end),yy(1:end-2)'-yy(3:end)'];    
0171     dangle = ((d1(1:end-2,1).*d1(3:end,2)...
0172         -d1(1:end-2,2).*d1(3:end,1))./...
0173         (norm2(d1(1:end-2,:)).*norm2(d1(3:end,:))));
0174     
0175     stm1 = max( dangle.^2 );    
0176     angle = cumsum(dangle)';
0177     na = numel(angle)-1;
0178     
0179     s = (0:na)/na*pi/2;
0180     
0181     try
0182         sa1 = mean( angle.*sin(s  ) );
0183         sa2 = mean( angle.*sin(s*2) );
0184         sa3 = mean( angle.*sin(s*3) );
0185         sa4 = mean( angle.*sin(s*4) );
0186     catch ME
0187         printError(ME)
0188     end
0189     
0190     
0191     if debug
0192         ss_ = size(d1);
0193         for ii = 1:ss_(1)
0194             plot( xx(ii)+[0,-d1(ii,1)], yy(ii)+[0,-d1(ii,2)], 'c.' );
0195         end
0196     end
0197     
0198     
0199     % Calculate the end score
0200     yy = [(ymin+dr),(ymax-dr)];
0201     xx = widthWindow(round(yy));
0202     r = ((xx)/2);
0203     r(r<1)=1;
0204     
0205     yyd = [(ymin+r(1)+[-1,0,1]),(ymax-r(2)+[-1,0,1])];
0206     yy  = round(yyd);
0207     
0208     [X] = meshgrid(x,yy);
0209     
0210     tmp = imRot(yy,:);
0211     xx = sum( X.*tmp,2 )./sum(tmp,2);
0212     xx = [mean(xx(1:3)),mean(xx(4:6))];
0213     yyd = yyd([2,5]);
0214     
0215     
0216     if debug
0217         phi = 0:.1:2*pi;
0218         
0219         plot(xx,yyd,'go');
0220         xxx = r(1)*cos(phi);
0221         yyy = r(1)*sin(phi);
0222         plot(xx(1)+xxx,yy(1)+yyy,'g:')
0223         
0224         xxx = r(2)*cos(phi);
0225         yyy = r(2)*sin(phi);
0226         plot(xx(2)+xxx,yyd(2)+yyy,'g:')
0227         
0228         plot( [1,ss(2)],[ymax,ymax],':r')
0229         plot( [1,ss(2)],[ymin,ymin],':r')
0230         
0231     end
0232     
0233     xcen = xx;
0234     ycen = yyd;
0235     
0236     x0 = xcen(1);
0237     y0 = ycen(1);
0238     r0 = r(1);
0239     
0240     yy = max(1,round(y0-r0-2)):min(round(y0+r0+2),ss(1));
0241     xx = max(1,round(x0-r0-2)):min(round(x0+r0+2),ss(2));
0242     
0243     [X,Y] = meshgrid(xx,yy);
0244     tmp = imRot(yy,xx);
0245     
0246     R2 = (X-x0).^2 + (Y-y0).^2 - r0^2;
0247     
0248     ind_over = find( ((R2 > 0).*tmp).*(d1(1,1)*(X-x0) + d1(1,2)*(Y-y0) > 0));
0249     ind_under = find( ((R2 < 0).*(tmp<1)).*(d1(1,1)*(X-x0) + d1(1,2)*(Y-y0) > 0));
0250     
0251     
0252     if debug
0253         plot(X(ind_over),Y(ind_over),'r.');
0254         plot(X(ind_under),Y(ind_under),'b.');
0255     end
0256     
0257     RoundIndOverTop  = sum( tmp(ind_over).*R2(ind_over))/r0^4;
0258     RoundIndUnderTop = sum( (1-tmp(ind_under)).*abs(R2(ind_under)))/r0^4;
0259     
0260     x0 = xcen(2);
0261     y0 = ycen(2);
0262     r0 = r(2);
0263     
0264     yy = max(1,round(y0-r0-2)):min(round(y0+r0+2),ss(1));
0265     xx = max(1,round(x0-r0-2)):min(round(x0+r0+2),ss(2));
0266     
0267     [X,Y] = meshgrid(xx,yy);
0268     tmp = imRot(yy,xx);
0269     
0270     R2 = (X-x0).^2 + (Y-y0).^2 - r0^2;
0271     
0272     ind_over = find( ((R2 > 0).*tmp).*(d1(1,1)*(X-x0) + d1(1,2)*(Y-y0) < 0));
0273     ind_under = find( ((R2 < 0).*(tmp<1)).*(d1(1,1)*(X-x0) + d1(1,2)*(Y-y0) < 0));
0274     
0275     if debug
0276         plot(X(ind_over),Y(ind_over),'r.');
0277         plot(X(ind_under),Y(ind_under),'b.');
0278     end
0279     
0280     
0281     RoundIndOverBot  = sum( tmp(ind_over).*R2(ind_over))/r0^4;
0282     RoundIndUnderBot = sum( (1-tmp(ind_under)).*abs(R2(ind_under)))/r0^4;
0283     
0284 else
0285     %'fail'
0286     stm1 = 0;
0287     sa1 = 0;
0288     sa2 = 0;
0289     sa3 = 0;
0290     sa4 = 0;
0291     RoundIndOverTop  = 0;
0292     RoundIndUnderTop = 0;
0293     RoundIndOverBot  = 0;
0294     RoundIndUnderBot = 0;
0295 end
0296 
0297 RoundIndOver  = RoundIndOverTop  + RoundIndOverBot;
0298 RoundIndUnder = RoundIndUnderTop + RoundIndUnderBot;
0299 
0300 if L1 < 1
0301     L1 = 1;
0302 end
0303 
0304 
0305 info = [L1, ...
0306     L2mean, ...
0307     Lneck, ...
0308     L2max, ...
0309     L2v, ...
0310     s1, ...
0311     s2, ...
0312     s3, ...
0313     RoundIndOver, ...
0314     RoundIndUnder, ...
0315     props.Area, ....
0316     stm1, ...
0317     sa1, ...
0318     sa2, ...
0319     sa3, ...
0320     sa4, ...
0321     1/L1, ...
0322     1/L2mean,...
0323     minWidthEnd,...
0324     sqmin,...
0325     sqmax];
0326 
0327 info(isnan(info)) = 0;
0328 
0329 if debug
0330     drawnow;
0331     keyboard;
0332 end
0333 
0334 end
0335 
0336 function tmp = norm2( dd )
0337 
0338 tmp = sqrt(sum(dd.*dd,2));
0339 
0340 end