Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse code

Integrated min peak intensity into dotfinder, changed z-stack method …

…for CheckDotUpDown. Fixed bug in lsmreader which created offset errors
  • Loading branch information...
commit 8400ff2b7f338fc798db777e2a1da81c4c9554a5 1 parent e17fee1
AlistairMobile authored November 22, 2011
120  CheckDotUpDown.m
@@ -24,10 +24,6 @@
24 24
 % use linear indexing of all dots
25 25
 tic 
26 26
 
27  
-     
28  
-% intype = 'uint8';
29  
-
30  
-disp('connecting dots in Z...') 
31 27
 %  plotdata = 1;  h = hs; w = ws; getpreciseZ =1; ovlap = 2; intype = 'uint16';  
32 28
 % DotData = DotData1; Inds = Inds1; Ints = Ints1; DotLabels = DotLabels1; 
33 29
 % DotData1 = DotData; Inds1 = Inds; Ints1 = Ints; DotLabels1 = DotLabels; 
@@ -48,7 +44,7 @@
48 44
     dotzpos{z} = z*ones(dotsinlayer(z),1);
49 45
 end
50 46
 dotC = [cell2mat(DotData'), cell2mat(dotzpos)];  
51  
- % clear DotData; 
  47
+clear DotData; 
52 48
 
53 49
 
54 50
 NDots = length(dotC); % total number of dots;
@@ -60,7 +56,7 @@
60 56
         % in 2 different data structures.  Could build this guy to start
61 57
         % with.  
62 58
 
63  
-       
  59
+%% Troubleshooting plots (3D dot scatter)       
64 60
 % figure(10); clf;  plot3(dotC(:,1),dotC(:,2),dotC(:,3),'w.');
65 61
 % % for n=1:NDots
66 62
 % %     z = dotC(n,3); 
@@ -88,7 +84,8 @@
88 84
 %  hold on; plot(DotData{21}(:,1),DotData{21}(:,2),'c.');
89 85
 % hold on; plot(DotData{22}(:,1),DotData{22}(:,2),'b.');
90 86
         
91  
-   
  87
+%%
  88
+disp('connecting dots in Z...') 
92 89
 DotConn = zeros(2*NDots,Zs,'single'); % empty connectivity matrix for all dots;  as a uint16 restricts this to 65,536 dots per image.  
93 90
 ConnInt = zeros(2*NDots,Zs,intype); 
94 91
 LayerJoin = false(2*NDots,Zs); 
@@ -122,8 +119,6 @@
122 119
 %             Iin_z = imreadfast(im_folder{Z});  
123 120
 %             I_test = Iin_z(xp1:xp2,yp1:yp2,mRNAchn+1);  imagesc(I_test); colorbar; colormap hot;
124 121
 %              hold on; plot(dotC(indsT,1),dotC(indsT,2),'c+');
125  
-            
126  
-   
127 122
          
128 123
          DotConn(2*st1_dot_num+1:2:2*(st1_dot_num + dotsinlayer(Z)),z) =  indsT; % STORE in DotConn matrix the indices 
129 124
          ConnInt(2*st1_dot_num+1:2:2*(st1_dot_num + dotsinlayer(Z)),z) = Ints{z}(inds_zin1+1); % Ints{z}(inds_zin1+1); ; %  % store actual intensities.  
@@ -154,18 +149,6 @@
154 149
 %     liny = dotC(nonzeros(DotConn(n,:)),2) ;  
155 150
 %     plot(linx,liny,'w');
156 151
 % end
157  
-
158  
-%    for z=1:Zs
159  
-%             text(DotData2{z}(:,1),DotData2{z}(:,2),[num2str(z)],'color','w','FontSize',8);
160  
-%            % text(DotData1{z}(:,1),DotData1{z}(:,2),[num2str(z)],'color','m','FontSize',8);
161  
-%    end
162  
-%         
163  
-
164  
-% % Flag specific points of interest.  
165  
-% p1 = (1089+1)/2; p2 = (525);
166  
-% figure(2); hold on; plot(dotC(p1,1),dotC(p1,2),'mo','MarkerSize',15); 
167  
-% plot(dotC(p2,1),dotC(p2,2),'m*','MarkerSize',15);
168  
-   
169 152
       
170 153
    
171 154
 
@@ -232,8 +215,6 @@
232 215
         labeled = bwlabel(mask);
233 216
         R = regionprops(labeled,ConnInt_T,'WeightedCentroid');
234 217
         tcent =  reshape([R.WeightedCentroid],2,length(R))';
235  
-
236  
-      %  cent = [cent; tcent(:,1),(k-1)*stp + tcent(:,2)];
237 218
         Cent{k} = [tcent(:,1),(k-1)*stp + tcent(:,2)];
238 219
         
239 220
         if plotdata == 1;
@@ -246,8 +227,8 @@
246 227
             title('Cross-Section of all dots'); 
247 228
         end
248 229
     end
249  
-
250 230
    
  231
+    % Final output is just DotConn split up by mask
251 232
        masked_inds(1+(k-1)*stp:min(k*stp,2*NDots),:) =  mask.*DotConn(1+(k-1)*stp:min(k*stp,2*NDots),:)  ;
252 233
 %      
253 234
        
@@ -256,6 +237,8 @@
256 237
        end
257 238
 end
258 239
 
  240
+% figure(3); clf; imagesc(masked_inds); colormap hot;
  241
+
259 242
 if plotdata == 1
260 243
     cent = cell2mat(Cent); 
261 244
 end
@@ -266,42 +249,65 @@
266 249
 
267 250
 
268 251
   %% Troubleshooting
269  
-% tic
270  
-% 
271 252
 % figure(2); clf; 
272 253
 % imagesc(Imax_dots);   hold on;
273 254
 % plot(dotC(:,1),dotC(:,2),'c+');
274 255
 
275 256
 %%
  257
+
  258
+%%  Free up some Memory
  259
+%clear  ConnInt_T DotConn LayerJoin  mask  ConnInt
  260
+% 
  261
+
  262
+%% New method for removing duplicates
  263
+
  264
+% masked inds has continuous blocks of connected dots, such that
  265
+% dotC(masked_inds) gives the the x and y coordinates of of the N
  266
+% associated dots with any given mRNA spot.  
  267
+
  268
+% we first pulll out these separate clusters, then we loop through strings
  269
+% of dots, comparing them to all other strings of dots, to remove the ones
  270
+% which occur twice. 
  271
+
  272
+tic
  273
+disp('Building unique-spheres from cross-section disks...');
  274
+
276 275
 Linx = cell(NDots,1); Liny = cell(NDots,1); 
277 276
 for n=1:2:2*NDots  
278  
-     % text(dotC((n+1)/2,1)+1,dotC((n+1)/2,2),['   ', num2str( (n+1)/2 )],'color','c','FontSize',8);
279  
-    % text(dotC((n+1)/2,1)+1,dotC((n+1)/2,2),['   ', num2str( dotC((n+1)/2,3) )],'color','c','FontSize',8);
280 277
     Linx{(n+1)/2} = dotC(nonzeros(masked_inds(n,:)),1);
281 278
     Liny{(n+1)/2} = dotC(nonzeros(masked_inds(n,:)),2);
282 279
     % % troubleshooting
283  
-    %plot(Linx{(n+1)/2},Liny{(n+1)/2},'c');    
  280
+    %plot(Linx{(n+1)/2},Liny{(n+1)/2},'c'); 
  281
+    % text(dotC((n+1)/2,1)+1,dotC((n+1)/2,2),['   ', num2str(dotC((n+1)/2,3) )],'color','c','FontSize',8);
284 282
 end
  283
+%%
285 284
 
286  
-
287  
-% dup = cell(length(Linx),1);
288 285
 for j=1:length(Linx); % j =1
  286
+         if isempty(Linx{j})==1;
  287
+             continue; 
  288
+         end
289 289
     for k=1:length(Linx) % k =3
290  
-        if j~=k
291  
-           if isempty(Linx{k})==0 && isempty(Linx{j})==0 
292  
-               if length(intersect(Linx{k}, Linx{j})) > 3
293  
-                 %  dup{j} = [dup{j},k];
  290
+         if isempty(Linx{k})==1;
  291
+             continue; 
  292
+         end
  293
+            if j~=k
  294
+               if (Linx{k}(1) - Linx{j}(1) + Liny{k}(1) - Liny{j}(1))^2 <  .01
  295
+                   % length(intersect(Linx{k}, Linx{j})) > 1  % this is the slow step
294 296
                   Linx{k} = [];
295  
-                 %  disp('duplicate removed');
  297
+                  % disp('duplicate removed'); 
296 298
                end
297  
-           end
298  
-        end
  299
+            end
299 300
     end
300  
- %   unique_ind(j)median(dup{j})
301 301
 end
302  
-old_dotC = dotC;
303  
-dotC = old_dotC; 
304  
-%%
  302
+
  303
+%  Convert string of dots to x-y-z of center dot. 
  304
+
  305
+% The median x,y position is taken as the true center since true 3D Gaussian
  306
+% fitting is too slow.
  307
+% the z position is either determined from the precise center of mass
  308
+% fitting of the trace if get-preciseZ is active, otherwise it is chosen as
  309
+% the middle of the stack (position where dot is first detected + half).  
  310
+
305 311
 unique_inds = find(~cellfun('isempty',Linx));
306 312
 unique_dotX = [Linx(~cellfun('isempty',Linx))];
307 313
 unique_dotY = [Liny( ~cellfun('isempty',Linx))];
@@ -309,9 +315,9 @@
309 315
 Ndots = length(unique_dotX);
310 316
 New_dotC = zeros(Ndots,3); 
311 317
 for k =1:Ndots
312  
-    phalf = round(Ndots/2);
313  
-    New_dotC(k,1) = unique_dotX{k}(phalf);
314  
-    New_dotC(k,2) = unique_dotY{k}(phalf);
  318
+    phalf = round(length(unique_dotX{k}/2));
  319
+    New_dotC(k,1) = median(unique_dotX{k}); % unique_dotX{k}(phalf);
  320
+    New_dotC(k,2) = median(unique_dotY{k}); % unique_dotY{k}(phalf);
315 321
     
316 322
     if plotdata == 1 && getpreciseZ == 1
317 323
         New_dotC(k,3) = cent(cent(:,2) == 2*unique_inds(k)-1,1);
@@ -320,23 +326,17 @@
320 326
     end
321 327
 end
322 328
 
  329
+disp([num2str(Ndots),' total spheres found']);
  330
+
  331
+toc
  332
+
323 333
 % 
324 334
 % figure(2); clf; 
325 335
 % imagesc(Imax_dots);   hold on;
326  
-% plot(New_dotC(:,1),New_dotC(:,2),'c+');
  336
+% plot(dotC(:,1),dotC(:,2),'c+');
  337
+% plot(New_dotC(:,1),New_dotC(:,2),'w.','MarkerSize',30);
327 338
 
328 339
 
329  
-  
330  
-% 
331  
-% hold on;
332  
-%    for z=1:Zs
333  
-%             text(DotData2{z}(:,1),DotData2{z}(:,2),[num2str(z)],'color','w','FontSize',8);
334  
-%            % text(DotData1{z}(:,1),DotData1{z}(:,2),[num2str(z)],'color','m','FontSize',8);
335  
-%    end
336  
-%         
337  
-% 
338  
-% 
339  
-% toc
340 340
 
341 341
 
342 342
 
@@ -346,10 +346,8 @@
346 346
 
347 347
 
348 348
 
349  
-%%
350  
-%clear  ConnInt_T DotConn LayerJoin  mask  ConnInt
351  
-% 
352  
-% %%  Remove duplicates
  349
+
  350
+%% Old Method to Remove duplicates
353 351
 % 
354 352
 % tic
355 353
 % disp('Building spheres from cross-section disks...');
183  CheckDotUpDown.m~
@@ -10,7 +10,7 @@
10 10
 % if getpreciseZ is off, the first layer in which the dot occurs is used as
11 11
 % the z postion, rather than the brightest layer  
12 12
 
13  
-function dotC = CheckDotUpDown(DotLabels,DotData,Inds,Ints,plotdata,getpreciseZ,consec_layers,ovlap,xp1,xp2,yp1,yp2,intype)
  13
+function New_dotC = CheckDotUpDown(DotLabels,DotData,Inds,Ints,plotdata,getpreciseZ,consec_layers,ovlap,xp1,xp2,yp1,yp2,intype)
14 14
 
15 15
 %% Updates
16 16
 % Rewritten 03/07/11 to convert more things to uint16 / uint8 to save
@@ -24,15 +24,12 @@ function dotC = CheckDotUpDown(DotLabels,DotData,Inds,Ints,plotdata,getpreciseZ,
24 24
 % use linear indexing of all dots
25 25
 tic 
26 26
 
27  
-     
28  
-% intype = 'uint8';
29  
-
30  
-disp('connecting dots in Z...') 
31 27
 %  plotdata = 1;  h = hs; w = ws; getpreciseZ =1; ovlap = 2; intype = 'uint16';  
32 28
 % DotData = DotData1; Inds = Inds1; Ints = Ints1; DotLabels = DotLabels1; 
33 29
 % DotData1 = DotData; Inds1 = Inds; Ints1 = Ints; DotLabels1 = DotLabels; 
34 30
  % DotData = DotData(mRNAchn,:); Inds = Inds(mRNAchn,:);  Ints = Ints(mRNAchn,:); DotLabels = DotLabels(mRNAchn,:);
35  
-  
  31
+ % mRNAchn = 2; 
  32
+ 
36 33
  comp_onVoff = 0;
37 34
 hs = yp2 - yp1+1; 
38 35
 ws = xp2 - xp1+1;
@@ -47,7 +44,7 @@ for z = 1:Zs
47 44
     dotzpos{z} = z*ones(dotsinlayer(z),1);
48 45
 end
49 46
 dotC = [cell2mat(DotData'), cell2mat(dotzpos)];  
50  
-% clear DotData; 
  47
+clear DotData; 
51 48
 
52 49
 
53 50
 NDots = length(dotC); % total number of dots;
@@ -59,7 +56,7 @@ disp(['Total dots = ',num2str(NDots)]);
59 56
         % in 2 different data structures.  Could build this guy to start
60 57
         % with.  
61 58
 
62  
-       
  59
+%% Troubleshooting plots (3D dot scatter)       
63 60
 % figure(10); clf;  plot3(dotC(:,1),dotC(:,2),dotC(:,3),'w.');
64 61
 % % for n=1:NDots
65 62
 % %     z = dotC(n,3); 
@@ -87,7 +84,8 @@ disp(['Total dots = ',num2str(NDots)]);
87 84
 %  hold on; plot(DotData{21}(:,1),DotData{21}(:,2),'c.');
88 85
 % hold on; plot(DotData{22}(:,1),DotData{22}(:,2),'b.');
89 86
         
90  
-   
  87
+%%
  88
+disp('connecting dots in Z...') 
91 89
 DotConn = zeros(2*NDots,Zs,'single'); % empty connectivity matrix for all dots;  as a uint16 restricts this to 65,536 dots per image.  
92 90
 ConnInt = zeros(2*NDots,Zs,intype); 
93 91
 LayerJoin = false(2*NDots,Zs); 
@@ -121,8 +119,6 @@ for Z = 1:Zs % The primary layer Z = 8
121 119
 %             Iin_z = imreadfast(im_folder{Z});  
122 120
 %             I_test = Iin_z(xp1:xp2,yp1:yp2,mRNAchn+1);  imagesc(I_test); colorbar; colormap hot;
123 121
 %              hold on; plot(dotC(indsT,1),dotC(indsT,2),'c+');
124  
-            
125  
-   
126 122
          
127 123
          DotConn(2*st1_dot_num+1:2:2*(st1_dot_num + dotsinlayer(Z)),z) =  indsT; % STORE in DotConn matrix the indices 
128 124
          ConnInt(2*st1_dot_num+1:2:2*(st1_dot_num + dotsinlayer(Z)),z) = Ints{z}(inds_zin1+1); % Ints{z}(inds_zin1+1); ; %  % store actual intensities.  
@@ -140,8 +136,8 @@ toc
140 136
 
141 137
 
142 138
 %%  Trouble-shooting: Draw lines between connected dots
143  
-
144  
-
  139
+% 
  140
+% 
145 141
 figure(1); clf; 
146 142
   Imax = imread([rawfolder,stackfolder,'max_',fname,'_',emb,'.tif']); 
147 143
             Imax_dots = 3*Imax(xp1:xp2,yp1:yp2,1:3);  
@@ -153,18 +149,6 @@ for n=1:2*NDots
153 149
     liny = dotC(nonzeros(DotConn(n,:)),2) ;  
154 150
     plot(linx,liny,'w');
155 151
 end
156  
-
157  
-%    for z=1:Zs
158  
-%             text(DotData2{z}(:,1),DotData2{z}(:,2),[num2str(z)],'color','w','FontSize',8);
159  
-%            % text(DotData1{z}(:,1),DotData1{z}(:,2),[num2str(z)],'color','m','FontSize',8);
160  
-%    end
161  
-%         
162  
-
163  
-% % Flag specific points of interest.  
164  
-% p1 = (1089+1)/2; p2 = (525);
165  
-% figure(2); hold on; plot(dotC(p1,1),dotC(p1,2),'mo','MarkerSize',15); 
166  
-% plot(dotC(p2,1),dotC(p2,2),'m*','MarkerSize',15);
167  
-   
168 152
       
169 153
    
170 154
 
@@ -231,8 +215,6 @@ for k=1:Nsects
231 215
         labeled = bwlabel(mask);
232 216
         R = regionprops(labeled,ConnInt_T,'WeightedCentroid');
233 217
         tcent =  reshape([R.WeightedCentroid],2,length(R))';
234  
-
235  
-      %  cent = [cent; tcent(:,1),(k-1)*stp + tcent(:,2)];
236 218
         Cent{k} = [tcent(:,1),(k-1)*stp + tcent(:,2)];
237 219
         
238 220
         if plotdata == 1;
@@ -245,8 +227,8 @@ for k=1:Nsects
245 227
             title('Cross-Section of all dots'); 
246 228
         end
247 229
     end
248  
-
249 230
    
  231
+    % Final output is just DotConn split up by mask
250 232
        masked_inds(1+(k-1)*stp:min(k*stp,2*NDots),:) =  mask.*DotConn(1+(k-1)*stp:min(k*stp,2*NDots),:)  ;
251 233
 %      
252 234
        
@@ -255,6 +237,8 @@ for k=1:Nsects
255 237
        end
256 238
 end
257 239
 
  240
+% figure(3); clf; imagesc(masked_inds); colormap hot;
  241
+
258 242
 if plotdata == 1
259 243
     cent = cell2mat(Cent); 
260 244
 end
@@ -265,81 +249,94 @@ toc
265 249
 
266 250
 
267 251
   %% Troubleshooting
268  
-% tic
  252
+% figure(2); clf; 
  253
+% imagesc(Imax_dots);   hold on;
  254
+% plot(dotC(:,1),dotC(:,2),'c+');
  255
+
  256
+%%
  257
+
  258
+%%  Free up some Memory
  259
+%clear  ConnInt_T DotConn LayerJoin  mask  ConnInt
269 260
 % 
270  
-figure(2); clf; 
271  
-imagesc(Imax_dots);   hold on;
272  
-plot(dotC(:,1),dotC(:,2),'c+');
273 261
 
  262
+%% New method for removing duplicates
274 263
 
275  
-Linx = cell(NDots,1); Liny = cell(NDots,1); 
  264
+% masked inds has continuous blocks of connected dots, such that
  265
+% dotC(masked_inds) gives the the x and y coordinates of of the N
  266
+% associated dots with any given mRNA spot.  
  267
+
  268
+% we first pulll out these separate clusters, then we loop through strings
  269
+% of dots, comparing them to all other strings of dots, to remove the ones
  270
+% which occur twice. 
  271
+
  272
+tic
  273
+disp('Building unique-spheres from cross-section disks...');
276 274
 
  275
+Linx = cell(NDots,1); Liny = cell(NDots,1); 
277 276
 for n=1:2:2*NDots  
278  
-     % text(dotC((n+1)/2,1)+1,dotC((n+1)/2,2),['   ', num2str( (n+1)/2 )],'color','c','FontSize',8);
279  
-     text(dotC((n+1)/2,1)+1,dotC((n+1)/2,2),['   ', num2str( dotC((n+1)/2,3) )],'color','c','FontSize',8);
280 277
     Linx{(n+1)/2} = dotC(nonzeros(masked_inds(n,:)),1);
281 278
     Liny{(n+1)/2} = dotC(nonzeros(masked_inds(n,:)),2);
282 279
     % % troubleshooting
283  
-    %plot(Linx{(n+1)/2},Liny{(n+1)/2},'c');    
  280
+    %plot(Linx{(n+1)/2},Liny{(n+1)/2},'c'); 
  281
+    % text(dotC((n+1)/2,1)+1,dotC((n+1)/2,2),['   ', num2str(dotC((n+1)/2,3) )],'color','c','FontSize',8);
284 282
 end
285  
-
286 283
 %%
287  
-% dup = cell(length(Linx),1);
  284
+
288 285
 for j=1:length(Linx); % j =1
  286
+         if isempty(Linx{j})==1;
  287
+             continue; 
  288
+         end
289 289
     for k=1:length(Linx) % k =3
290  
-        if j~=k
291  
-           if isempty(Linx{k})==0 && isempty(Linx{j})==0 
292  
-               if length(intersect(Linx{k}, Linx{j})) > 3
293  
-                 %  dup{j} = [dup{j},k];
  290
+         if isempty(Linx{k})==1;
  291
+             continue; 
  292
+         end
  293
+            if j~=k
  294
+               if (Linx{k}(1) - Linx{j}(1) + Liny{k}(1) - Liny{j}(1))^2 <  1
  295
+                   % length(intersect(Linx{k}, Linx{j})) > 1  % this is the slow step
294 296
                   Linx{k} = [];
295  
-                 %  disp('duplicate removed');
  297
+                  % disp('duplicate removed'); 
296 298
                end
297  
-           end
298  
-        end
  299
+            end
299 300
     end
300  
- %   unique_ind(j)median(dup{j})
301 301
 end
302  
-old_dotC = dotC;
303  
-dotC = old_dotC; 
304  
-%%
  302
+
  303
+%  Convert string of dots to x-y-z of center dot. 
  304
+
  305
+% The median x,y position is taken as the true center since true 3D Gaussian
  306
+% fitting is too slow.
  307
+% the z position is either determined from the precise center of mass
  308
+% fitting of the trace if get-preciseZ is active, otherwise it is chosen as
  309
+% the middle of the stack (position where dot is first detected + half).  
  310
+
305 311
 unique_inds = find(~cellfun('isempty',Linx));
306  
-unique_dotX = [Linx(~cellfun('isempty',Linx))]
307  
-unique_dotY = [Liny( ~cellfun('isempty',Linx))]
  312
+unique_dotX = [Linx(~cellfun('isempty',Linx))];
  313
+unique_dotY = [Liny( ~cellfun('isempty',Linx))];
308 314
 
309 315
 Ndots = length(unique_dotX);
310  
-
  316
+New_dotC = zeros(Ndots,3); 
311 317
 for k =1:Ndots
312  
-    phalf = round(Ndots/2);
313  
-    new_dotC(k,1) = unique_dotX{k}(phalf);
314  
-    new_dotC(k,2) = unique_dotY{k}(phalf);
  318
+    phalf = round(length(unique_dotX{k}/2));
  319
+    New_dotC(k,1) = median(unique_dotX{k}); % unique_dotX{k}(phalf);
  320
+    New_dotC(k,2) = median(unique_dotY{k}); % unique_dotY{k}(phalf);
315 321
     
316  
-    if getpreciseZ == 1
317  
-        zF = cent(cent(:,2) == unique_inds(k),1);
  322
+    if plotdata == 1 && getpreciseZ == 1
  323
+        New_dotC(k,3) = cent(cent(:,2) == 2*unique_inds(k)-1,1);
318 324
     else
319  
-        zF = dotC(unique_inds(k),3) + phalf;
  325
+        New_dotC(k,3) = dotC(unique_inds(k),3) + phalf;
320 326
     end
321  
-  
322  
-%cent(~cellfun('isempty',Linx),1)
323  
-xF = dotC(~cellfun('isempty',Linx),1)
324  
-yF = dotC(~cellfun('isempty',Linx),2)
325  
-dotC(~cellfun('isempty',Linx),3)
  327
+end
  328
+
  329
+disp([num2str(Ndots),' total spheres found']);
  330
+
  331
+toc
326 332
 
  333
+% 
327 334
 figure(2); clf; 
328 335
 imagesc(Imax_dots);   hold on;
329  
-plot(xF,yF,'c+');
  336
+plot(dotC(:,1),dotC(:,2),'c+');
  337
+plot(New_dotC(:,1),New_dotC(:,2),'w.','MarkerSize',30);
330 338
 
331 339
 
332  
-  
333  
-% 
334  
-% hold on;
335  
-%    for z=1:Zs
336  
-%             text(DotData2{z}(:,1),DotData2{z}(:,2),[num2str(z)],'color','w','FontSize',8);
337  
-%            % text(DotData1{z}(:,1),DotData1{z}(:,2),[num2str(z)],'color','m','FontSize',8);
338  
-%    end
339  
-%         
340  
-% 
341  
-% 
342  
-% toc
343 340
 
344 341
 
345 342
 
@@ -349,10 +346,8 @@ plot(xF,yF,'c+');
349 346
 
350 347
 
351 348
 
352  
-%%
353  
-%clear  ConnInt_T DotConn LayerJoin  mask  ConnInt
354  
-% 
355  
-% %%  Remove duplicates
  349
+
  350
+%% Old Method to Remove duplicates
356 351
 % 
357 352
 % tic
358 353
 % disp('Building spheres from cross-section disks...');
@@ -433,21 +428,21 @@ plot(xF,yF,'c+');
433 428
  
434 429
     %%      
435 430
           
436  
-%New_dotC = dotC(~remove_dot,:);
437  
-dotC = dotC(~remove_dot,:);
438  
-[N_dots, jnk] = size(dotC); 
439  
-%N_dots = NDots - sum(remove_dot) % sum(stacked_dots)
440  
-%N_dots = length(New_dotC); 
441  
-disp(['Counted ',num2str(N_dots),' spheres']); 
442  
-
443  
-
444  
-
445  
-if comp_onVoff == 1
446  
-    figure(7); clf;  subplot(2,1,1);
447  
-    hist(log(nonzeros(masked_ints(~remove_dot,:)))); title('intensities of kept dots');
448  
-    subplot(2,1,2);
449  
-    hist(log(nonzeros(masked_ints(logical(remove_dot),:)))); title('intensities of removed dots');
450  
-end
  431
+% %New_dotC = dotC(~remove_dot,:);
  432
+% dotC = dotC(~remove_dot,:);
  433
+% [N_dots, jnk] = size(dotC); 
  434
+% %N_dots = NDots - sum(remove_dot) % sum(stacked_dots)
  435
+% %N_dots = length(New_dotC); 
  436
+% disp(['Counted ',num2str(N_dots),' spheres']); 
  437
+% 
  438
+% 
  439
+% 
  440
+% if comp_onVoff == 1
  441
+%     figure(7); clf;  subplot(2,1,1);
  442
+%     hist(log(nonzeros(masked_ints(~remove_dot,:)))); title('intensities of kept dots');
  443
+%     subplot(2,1,2);
  444
+%     hist(log(nonzeros(masked_ints(logical(remove_dot),:)))); title('intensities of removed dots');
  445
+% end
451 446
 
452 447
 
453 448
 
53  Unsupervised_DotFinding.m
@@ -8,12 +8,12 @@
8 8
 
9 9
 tot_time = tic;
10 10
 % Input options 
11  
-old_lab = 0;  Es = 0;  ver = ''; %   '_v4'; % '_v3';% '_v2';
  11
+old_lab = 0;  Es = 0;  ver = '';
12 12
 folder = '/home/alistair/Documents/Research/Projects/mRNA_counting/Data/2011-05-22/'; %2011-06-20/'; %  2011-04_and_earlier/'; % % 2011-05-22/'; % 2011-06-20/'; %   '/Users/alistair/Documents/Berkeley/Levine_Lab/Projects/Enhancer_Modeling/Data/'; 
13 13
 rawfolder = '/home/alistair/Documents/Research/Raw_Data_Temp/2011-05-22/'; %2011-06-20/'; %  2011-04_and_earlier/'; %'; % 2011-06-20/'; %  '/Volumes/Data/Lab Data/Raw_Data/02-17-11/'; %%   %
14  
-
15  
-stackfolder = 's07_MP08/'; % 's07_MP05Hz/';% 's04_MP10/';%   'MP07Hz/'; %  'MP12Hz/'; %    's05_MP06/';% 's02_MP01/';% 's01_MP09/';%   'sna2.8Hz/' ;%'s06_MP10_sna18/'; %'s21_MP07/';% 'MP07Hz/';% 's11_G4B/' %  's06_MP10_sna18/'; % %'s10_bcd1x/';%  's11_bcd6x/'; %'s14_comp_cntrl/'; % 's12_cntrl_2label/'; %'MP02_22C/'; %'MP01_22C/'; % 'MGa1x/'; % 'MP10_22C/'; %'MP05_22C/'; %'YW_ths_sog/'; % 'MP10_22C/'; %  % 'MP09_22C/'; % 'MGa2x/'; % 'MGa1x/'; % 'MGa2x/'; % 'MP10_22C_sna_y_c/'; %
16  
-fname =  'MP08Hz_snaD_22C_b'; % 's07_MP05Hz_22C'; ver = '_v2'; % 'MP10Hz_c'; %'MP07Hz_snaD_22C_b' ; ver = '_v3';%  'MP12Hz_snaD_22C_b'; %    's04_MP10Hz'; % 's05_MP06Hz_b';% 's02_MP01_Hz_22C_b'; % 's01_MP09_Hz_22C_c'; %'sna2.8Hz_snaD_22C'; % 's06_MP10_sna18_b'; % 'MP07het_snaD_22C'; %  'MP07Hz_snaD_22C';%'s11_G4B_LacZ';% 's06_MP10_sna18_b'; % 's05_MP06Hz'; %   %'s10_bcd1x';% 's11_bcd6x'; % 's14_comp_cntrl'; Es =1; % 's12_cntrl_2label'; Es = 1; % 'MP09_22C_hb_y_f'; Es = 7; %  'MP02_22C_hb_y'; Es = 9; % 'MP02_22C_hb_y_b'; Es = 10; %  % 'MP01_22C_hb_y_f'; Es = 12; % 'MP01_22C_hb_y_c'; Es = 10; % 'MP01_22C_hb_y'; Es = 13; % 'MGa1x_LacZ_b'; Es = 12; %  'MP10_22C_sna_y_e'; Es = 12; %  'MP05_22C_sna_y_c'; Es =7; %  'MP10_22C_sna_y_d3'; Es = 1;  %'YW_ths_sog'; Es = 12;  % % 'MP09_22C_hb_y_e'; Es = 10; % 'MP09_22C_hb_y_d'; Es=11; % 'MGa2x_LacZ_sna_b'; Es = 10; % 'MP10_22C_sna_y_d';   % 'MGa_LacZ'; %'MGa2x_LacZ_sna'; %'MP10_22C_sna_y_c'; old_lab = 1;  % 'MP05_22C_sna_y'; old_lab = 1; % 
  14
+rawfolder =  '/media/GRAID/Raw_Data/2011-05-22/'
  15
+stackfolder =  's05_MP06/';% 's07_MP08/'; % 's07_MP05Hz/';% 's04_MP10/';%   'MP07Hz/'; %  'MP12Hz/'; %    's02_MP01/';% 's01_MP09/';%   'sna2.8Hz/' ;%'s06_MP10_sna18/'; %'s21_MP07/';% 'MP07Hz/';% 's11_G4B/' %  's06_MP10_sna18/'; % %'s10_bcd1x/';%  's11_bcd6x/'; %'s14_comp_cntrl/'; % 's12_cntrl_2label/'; %'MP02_22C/'; %'MP01_22C/'; % 'MGa1x/'; % 'MP10_22C/'; %'MP05_22C/'; %'YW_ths_sog/'; % 'MP10_22C/'; %  % 'MP09_22C/'; % 'MGa2x/'; % 'MGa1x/'; % 'MGa2x/'; % 'MP10_22C_sna_y_c/'; %
  16
+fname =  's05_MP06Hz_b', ver = '_v3'% 'MP08Hz_snaD_22C_b'; % 's07_MP05Hz_22C'; ver = '_v2'; % 'MP10Hz_c'; %'MP07Hz_snaD_22C_b' ; ver = '_v3';%  'MP12Hz_snaD_22C_b'; %    's04_MP10Hz'; % 's02_MP01_Hz_22C_b'; % 's01_MP09_Hz_22C_c'; %'sna2.8Hz_snaD_22C'; % 's06_MP10_sna18_b'; % 'MP07het_snaD_22C'; %  'MP07Hz_snaD_22C';%'s11_G4B_LacZ';% 's06_MP10_sna18_b'; % 's05_MP06Hz'; %   %'s10_bcd1x';% 's11_bcd6x'; % 's14_comp_cntrl'; Es =1; % 's12_cntrl_2label'; Es = 1; % 'MP09_22C_hb_y_f'; Es = 7; %  'MP02_22C_hb_y'; Es = 9; % 'MP02_22C_hb_y_b'; Es = 10; %  % 'MP01_22C_hb_y_f'; Es = 12; % 'MP01_22C_hb_y_c'; Es = 10; % 'MP01_22C_hb_y'; Es = 13; % 'MGa1x_LacZ_b'; Es = 12; %  'MP10_22C_sna_y_e'; Es = 12; %  'MP05_22C_sna_y_c'; Es =7; %  'MP10_22C_sna_y_d3'; Es = 1;  %'YW_ths_sog'; Es = 12;  % % 'MP09_22C_hb_y_e'; Es = 10; % 'MP09_22C_hb_y_d'; Es=11; % 'MGa2x_LacZ_sna_b'; Es = 10; % 'MP10_22C_sna_y_d';   % 'MGa_LacZ'; %'MGa2x_LacZ_sna'; %'MP10_22C_sna_y_c'; old_lab = 1;  % 'MP05_22C_sna_y'; old_lab = 1; % 
17 17
 mRNA_channels = 2;% 2; %  3; %  1; % total mRNA channels
18 18
 
19 19
 sname =  fname; % 'MP07het_snaD_22C_1';% '_1'; % additional label on slide. 
@@ -36,18 +36,18 @@
36 36
     Es = length(fields(Datas)) - 3;   % Number of Stacks
37 37
 end
38 38
 % ------- Option: Focus on subset of image: ------------------- %
39  
-%      m = 1/2048;  % .9;%    .7; % .5; .7; %   1/2048; % 
40  
-%    xp1= floor(h/2*m)+1; xp2 = floor(h/2*(2-m))+1;  yp1 = floor(w/2*m)+1;  yp2 = floor(w/2*(2-m))+1;
41  
-%    hs = yp2-yp1+1;     ws = xp2-xp1+1;
  39
+     m =   1/2048;  %    .7; % .5; .7; %   1/2048; % 
  40
+   xp1= floor(h/2*m)+1; xp2 = floor(h/2*(2-m))+1;  yp1 = floor(w/2*m)+1;  yp2 = floor(w/2*(2-m))+1;
  41
+   hs = yp2-yp1+1;     ws = xp2-xp1+1;
42 42
 
43  
-ws = 2048; hs = 2048; xp1 = 1; yp1 = 1;
44  
-xp2 = xp1 + ws -1; yp2 = yp1 + hs - 1; 
  43
+% ws = 2048; hs = 2048; xp1 = 1; yp1 = 1;
  44
+% xp2 = xp1 + ws -1; yp2 = yp1 + hs - 1; 
45 45
 disp(['Coordinates:  ', num2str(xp1), ' : ', num2str(xp2), ',   ' num2str(yp1), ' : ', num2str(yp2) ] );
46 46
 % ------------------------------------------------------------- %     
47 47
     
48 48
 
49 49
 % -------------- Graphing and Display Options ------------------ %
50  
-   show_projected = 0; % show max-project with all dots and linked dots. 
  50
+   show_projected = 1; % show max-project with all dots and linked dots. 
51 51
    plotdata = 0; % CheckDotUpDown display parameter
52 52
    plotZdata = 0 ;% show z-map of data
53 53
    showhist = 1; % show histogram of mRNA counts per cell. 
@@ -62,14 +62,17 @@
62 62
    % dotfinder's parameters 
63 63
     sigmaE = 3;%  IMPORTANT
64 64
     sigmaI = 4; % IMPORTANT
65  
-  %  min_int  = 0.04;    %  5    ;% .05 % not necessary Fix at Zero
  65
+    min_int  = 0.04;    %  5    ;% .05 % not necessary Fix at Zero
66 66
     FiltSize = 30;% 
67 67
     min_size = 30;% 
  68
+    min_peak = 600; %
68 69
    
69 70
   % sphere finding parameters
70 71
    getpreciseZ = 0;
71  
-  % consec_layers = 3;
72  
-   ovlap = 4; 
  72
+   consec_layers = 3;
  73
+   ovlap = 2;  
  74
+   % large ovlap yields confusing dots and then watershed splits these up
  75
+   % in weird dot-distructive ways
73 76
 %---------------------------------%
74 77
 
75 78
 
@@ -81,7 +84,7 @@
81 84
     disp('Running DotFinder1'); 
82 85
   
83 86
 %%
84  
-for e= 12:Es
  87
+for e= 1:Es
85 88
 %%
86 89
 disp('loading data...');
87 90
     tic 
@@ -91,7 +94,7 @@
91 94
     else
92 95
         emb = num2str(e);
93 96
     end
94  
-    disp(['analyzing embryo, ',emb,'...']);
  97
+  
95 98
    
96 99
     
97 100
     try load([rawfolder,stackfolder,fname,'_',emb,'_nucdata.mat']);
@@ -114,21 +117,11 @@
114 117
     end
115 118
    
116 119
     toc
117  
-    
  120
+      disp(['analyzing embryo, ',emb,'...']);
118 121
     
119 122
     for mRNAchn = 1:mRNA_channels % mRNAchn =2
120 123
         
121  
-         if mRNAchn == 1;
122  
-               min_int  = 0.05;  % just for speed 
123  
-               consec_layers = 3;
124  
-         elseif mRNAchn == 2
125  
-               min_int  = 0.02; %
126  
-                consec_layers = 3;
127  
-         elseif mRNAchn == 3
128  
-               min_int = .01;
129  
-               consec_layers = 3;
130  
-         end
131  
-        
  124
+           
132 125
             DotLabels= cell(1,Zs); 
133 126
             DotData = cell(1,Zs);    
134 127
             Inds = cell(1,Zs); 
@@ -140,7 +133,7 @@
140 133
                  im_folder{z} = [rawfolder,stackfolder,fname,'_',emb,'_z',num2str(z),'.tif'];
141 134
                  try
142 135
                  Iin_z = imreadfast(im_folder{z});       
143  
-                 [DotLabels{z},DotData{z},Inds{z},Ints{z}]  = dotfinder(Iin_z(xp1:xp2,yp1:yp2,mRNAchn),Ex,Ix,min_int,min_size);
  136
+                 [DotLabels{z},DotData{z},Inds{z},Ints{z}]  = dotfinder(Iin_z(xp1:xp2,yp1:yp2,mRNAchn),Ex,Ix,min_int,min_size,min_peak);
144 137
                  catch err
145 138
                      disp(err.message); 
146 139
                      Zs = z-1; 
@@ -174,7 +167,7 @@
174 167
             
175 168
             Imax_dots = Imax(xp1:xp2,yp1:yp2,mRNAchn);  
176 169
             figure(2); 
177  
-            Iout = figure(2);  clf;  imagesc(Imax_dots);
  170
+            Iout = figure(2);  clf;  imagesc(Imax_dots); colorbar;
178 171
             colordef black; set(gcf,'color','k'); 
179 172
             colormap hot; hold on;
180 173
             plot(  dotC(:,1),dotC(:,2),'w+','MarkerSize',14 );
@@ -219,7 +212,7 @@
219 212
         mRNA_den = zeros(1,Nnucs);  % store densities of mRNA per cell
220 213
         nuc_area = zeros(1,Nnucs); 
221 214
         if showim == 1
222  
-            Plot_mRNA = single(NucLabeled);
  215
+            Plot_mRNA = single(NucLabel);
223 216
         end
224 217
         for i=1:Nnucs; % i = 4
225 218
             nn = Nucs_list(i);
57  Unsupervised_DotFinding.m~
@@ -8,12 +8,12 @@ clear all;
8 8
 
9 9
 tot_time = tic;
10 10
 % Input options 
11  
-old_lab = 0;  Es = 0;  ver = ''; %   '_v4'; % '_v3';% '_v2';
12  
-folder = '/Users/alistair/Documents/Berkeley/Levine_Lab/Projects/mRNA_counting/Data/2011-05-22/'; %2011-06-20/'; %  2011-04_and_earlier/'; % % 2011-05-22/'; % 2011-06-20/'; %   '/Users/alistair/Documents/Berkeley/Levine_Lab/Projects/Enhancer_Modeling/Data/'; 
13  
-rawfolder = '/Volumes/Data/Lab Data/Raw_Data/2011-06-20/'; %  2011-04_and_earlier/'; %'; % 2011-06-20/'; %  '/Volumes/Data/Lab Data/Raw_Data/02-17-11/'; %%   %
14  
-
15  
-stackfolder = 's07_MP08/'; % 's07_MP05Hz/';% 's04_MP10/';%   'MP07Hz/'; %  'MP12Hz/'; %    's05_MP06/';% 's02_MP01/';% 's01_MP09/';%   'sna2.8Hz/' ;%'s06_MP10_sna18/'; %'s21_MP07/';% 'MP07Hz/';% 's11_G4B/' %  's06_MP10_sna18/'; % %'s10_bcd1x/';%  's11_bcd6x/'; %'s14_comp_cntrl/'; % 's12_cntrl_2label/'; %'MP02_22C/'; %'MP01_22C/'; % 'MGa1x/'; % 'MP10_22C/'; %'MP05_22C/'; %'YW_ths_sog/'; % 'MP10_22C/'; %  % 'MP09_22C/'; % 'MGa2x/'; % 'MGa1x/'; % 'MGa2x/'; % 'MP10_22C_sna_y_c/'; %
16  
-fname =  'MP08Hz_snaD_22C_b'; % 's07_MP05Hz_22C'; ver = '_v2'; % 'MP10Hz_c'; %'MP07Hz_snaD_22C_b' ; ver = '_v3';%  'MP12Hz_snaD_22C_b'; %    's04_MP10Hz'; % 's05_MP06Hz_b';% 's02_MP01_Hz_22C_b'; % 's01_MP09_Hz_22C_c'; %'sna2.8Hz_snaD_22C'; % 's06_MP10_sna18_b'; % 'MP07het_snaD_22C'; %  'MP07Hz_snaD_22C';%'s11_G4B_LacZ';% 's06_MP10_sna18_b'; % 's05_MP06Hz'; %   %'s10_bcd1x';% 's11_bcd6x'; % 's14_comp_cntrl'; Es =1; % 's12_cntrl_2label'; Es = 1; % 'MP09_22C_hb_y_f'; Es = 7; %  'MP02_22C_hb_y'; Es = 9; % 'MP02_22C_hb_y_b'; Es = 10; %  % 'MP01_22C_hb_y_f'; Es = 12; % 'MP01_22C_hb_y_c'; Es = 10; % 'MP01_22C_hb_y'; Es = 13; % 'MGa1x_LacZ_b'; Es = 12; %  'MP10_22C_sna_y_e'; Es = 12; %  'MP05_22C_sna_y_c'; Es =7; %  'MP10_22C_sna_y_d3'; Es = 1;  %'YW_ths_sog'; Es = 12;  % % 'MP09_22C_hb_y_e'; Es = 10; % 'MP09_22C_hb_y_d'; Es=11; % 'MGa2x_LacZ_sna_b'; Es = 10; % 'MP10_22C_sna_y_d';   % 'MGa_LacZ'; %'MGa2x_LacZ_sna'; %'MP10_22C_sna_y_c'; old_lab = 1;  % 'MP05_22C_sna_y'; old_lab = 1; % 
  11
+old_lab = 0;  Es = 0;  ver = '';
  12
+folder = '/home/alistair/Documents/Research/Projects/mRNA_counting/Data/2011-05-22/'; %2011-06-20/'; %  2011-04_and_earlier/'; % % 2011-05-22/'; % 2011-06-20/'; %   '/Users/alistair/Documents/Berkeley/Levine_Lab/Projects/Enhancer_Modeling/Data/'; 
  13
+rawfolder = '/home/alistair/Documents/Research/Raw_Data_Temp/2011-05-22/'; %2011-06-20/'; %  2011-04_and_earlier/'; %'; % 2011-06-20/'; %  '/Volumes/Data/Lab Data/Raw_Data/02-17-11/'; %%   %
  14
+rawfolder =  '/media/GRAID/Raw_Data/2011-05-22/'
  15
+stackfolder =  's05_MP06/';% 's07_MP08/'; % 's07_MP05Hz/';% 's04_MP10/';%   'MP07Hz/'; %  'MP12Hz/'; %    's02_MP01/';% 's01_MP09/';%   'sna2.8Hz/' ;%'s06_MP10_sna18/'; %'s21_MP07/';% 'MP07Hz/';% 's11_G4B/' %  's06_MP10_sna18/'; % %'s10_bcd1x/';%  's11_bcd6x/'; %'s14_comp_cntrl/'; % 's12_cntrl_2label/'; %'MP02_22C/'; %'MP01_22C/'; % 'MGa1x/'; % 'MP10_22C/'; %'MP05_22C/'; %'YW_ths_sog/'; % 'MP10_22C/'; %  % 'MP09_22C/'; % 'MGa2x/'; % 'MGa1x/'; % 'MGa2x/'; % 'MP10_22C_sna_y_c/'; %
  16
+fname =  's05_MP06Hz_b', ver = '_v3'% 'MP08Hz_snaD_22C_b'; % 's07_MP05Hz_22C'; ver = '_v2'; % 'MP10Hz_c'; %'MP07Hz_snaD_22C_b' ; ver = '_v3';%  'MP12Hz_snaD_22C_b'; %    's04_MP10Hz'; % 's02_MP01_Hz_22C_b'; % 's01_MP09_Hz_22C_c'; %'sna2.8Hz_snaD_22C'; % 's06_MP10_sna18_b'; % 'MP07het_snaD_22C'; %  'MP07Hz_snaD_22C';%'s11_G4B_LacZ';% 's06_MP10_sna18_b'; % 's05_MP06Hz'; %   %'s10_bcd1x';% 's11_bcd6x'; % 's14_comp_cntrl'; Es =1; % 's12_cntrl_2label'; Es = 1; % 'MP09_22C_hb_y_f'; Es = 7; %  'MP02_22C_hb_y'; Es = 9; % 'MP02_22C_hb_y_b'; Es = 10; %  % 'MP01_22C_hb_y_f'; Es = 12; % 'MP01_22C_hb_y_c'; Es = 10; % 'MP01_22C_hb_y'; Es = 13; % 'MGa1x_LacZ_b'; Es = 12; %  'MP10_22C_sna_y_e'; Es = 12; %  'MP05_22C_sna_y_c'; Es =7; %  'MP10_22C_sna_y_d3'; Es = 1;  %'YW_ths_sog'; Es = 12;  % % 'MP09_22C_hb_y_e'; Es = 10; % 'MP09_22C_hb_y_d'; Es=11; % 'MGa2x_LacZ_sna_b'; Es = 10; % 'MP10_22C_sna_y_d';   % 'MGa_LacZ'; %'MGa2x_LacZ_sna'; %'MP10_22C_sna_y_c'; old_lab = 1;  % 'MP05_22C_sna_y'; old_lab = 1; % 
17 17
 mRNA_channels = 2;% 2; %  3; %  1; % total mRNA channels
18 18
 
19 19
 sname =  fname; % 'MP07het_snaD_22C_1';% '_1'; % additional label on slide. 
@@ -36,18 +36,18 @@ if Es==0
36 36
     Es = length(fields(Datas)) - 3;   % Number of Stacks
37 37
 end
38 38
 % ------- Option: Focus on subset of image: ------------------- %
39  
-%      m = 1/2048;  % .9;%    .7; % .5; .7; %   1/2048; % 
40  
-%    xp1= floor(h/2*m)+1; xp2 = floor(h/2*(2-m))+1;  yp1 = floor(w/2*m)+1;  yp2 = floor(w/2*(2-m))+1;
41  
-%    hs = yp2-yp1+1;     ws = xp2-xp1+1;
  39
+     m = .9;%  1/2048;  %    .7; % .5; .7; %   1/2048; % 
  40
+   xp1= floor(h/2*m)+1; xp2 = floor(h/2*(2-m))+1;  yp1 = floor(w/2*m)+1;  yp2 = floor(w/2*(2-m))+1;
  41
+   hs = yp2-yp1+1;     ws = xp2-xp1+1;
42 42
 
43  
-ws = 2048; hs = 2048; xp1 = 1; yp1 = 1;
44  
-xp2 = xp1 + ws -1; yp2 = yp1 + hs - 1; 
  43
+% ws = 2048; hs = 2048; xp1 = 1; yp1 = 1;
  44
+% xp2 = xp1 + ws -1; yp2 = yp1 + hs - 1; 
45 45
 disp(['Coordinates:  ', num2str(xp1), ' : ', num2str(xp2), ',   ' num2str(yp1), ' : ', num2str(yp2) ] );
46 46
 % ------------------------------------------------------------- %     
47 47
     
48 48
 
49 49
 % -------------- Graphing and Display Options ------------------ %
50  
-   show_projected = 0; % show max-project with all dots and linked dots. 
  50
+   show_projected = 1; % show max-project with all dots and linked dots. 
51 51
    plotdata = 0; % CheckDotUpDown display parameter
52 52
    plotZdata = 0 ;% show z-map of data
53 53
    showhist = 1; % show histogram of mRNA counts per cell. 
@@ -62,14 +62,17 @@ disp(['Coordinates:  ', num2str(xp1), ' : ', num2str(xp2), ',   ' num2str(yp1),
62 62
    % dotfinder's parameters 
63 63
     sigmaE = 3;%  IMPORTANT
64 64
     sigmaI = 4; % IMPORTANT
65  
-  %  min_int  = 0.04;    %  5    ;% .05 % not necessary Fix at Zero
  65
+    min_int  = 0.04;    %  5    ;% .05 % not necessary Fix at Zero
66 66
     FiltSize = 30;% 
67 67
     min_size = 30;% 
  68
+    min_peak = 600; %
68 69
    
69 70
   % sphere finding parameters
70 71
    getpreciseZ = 0;
71  
-  % consec_layers = 3;
72  
-   ovlap = 4; 
  72
+   consec_layers = 3;
  73
+   ovlap = 2;  
  74
+   % large ovlap yields confusing dots and then watershed splits these up
  75
+   % in weird dot-distructive ways
73 76
 %---------------------------------%
74 77
 
75 78
 
@@ -81,7 +84,7 @@ disp(['Coordinates:  ', num2str(xp1), ' : ', num2str(xp2), ',   ' num2str(yp1),
81 84
     disp('Running DotFinder1'); 
82 85
   
83 86
 %%
84  
-for e= 12:Es
  87
+for e= 1:Es
85 88
 %%
86 89
 disp('loading data...');
87 90
     tic 
@@ -91,7 +94,7 @@ disp('loading data...');
91 94
     else
92 95
         emb = num2str(e);
93 96
     end
94  
-    disp(['analyzing embryo, ',emb,'...']);
  97
+  
95 98
    
96 99
     
97 100
     try load([rawfolder,stackfolder,fname,'_',emb,'_nucdata.mat']);
@@ -114,21 +117,11 @@ disp('loading data...');
114 117
     end
115 118
    
116 119
     toc
117  
-    
  120
+      disp(['analyzing embryo, ',emb,'...']);
118 121
     
119 122
     for mRNAchn = 1:mRNA_channels % mRNAchn =2
120 123
         
121  
-         if mRNAchn == 1;
122  
-               min_int  = 0.05;  % just for speed 
123  
-               consec_layers = 3;
124  
-         elseif mRNAchn == 2
125  
-               min_int  = 0.02; %
126  
-                consec_layers = 3;
127  
-         elseif mRNAchn == 3
128  
-               min_int = .01;
129  
-               consec_layers = 3;
130  
-         end
131  
-        
  124
+           
132 125
             DotLabels= cell(1,Zs); 
133 126
             DotData = cell(1,Zs);    
134 127
             Inds = cell(1,Zs); 
@@ -140,7 +133,7 @@ disp('loading data...');
140 133
                  im_folder{z} = [rawfolder,stackfolder,fname,'_',emb,'_z',num2str(z),'.tif'];
141 134
                  try
142 135
                  Iin_z = imreadfast(im_folder{z});       
143  
-                 [DotLabels{z},DotData{z},Inds{z},Ints{z}]  = dotfinder(Iin_z(xp1:xp2,yp1:yp2,mRNAchn),Ex,Ix,min_int,min_size);
  136
+                 [DotLabels{z},DotData{z},Inds{z},Ints{z}]  = dotfinder(Iin_z(xp1:xp2,yp1:yp2,mRNAchn),Ex,Ix,min_int,min_size,min_peak);
144 137
                  catch err
145 138
                      disp(err.message); 
146 139
                      Zs = z-1; 
@@ -174,7 +167,7 @@ disp('loading data...');
174 167
             
175 168
             Imax_dots = Imax(xp1:xp2,yp1:yp2,mRNAchn);  
176 169
             figure(2); 
177  
-            Iout = figure(2);  clf;  imagesc(Imax_dots);
  170
+            Iout = figure(2);  clf;  imagesc(Imax_dots); colorbar;
178 171
             colordef black; set(gcf,'color','k'); 
179 172
             colormap hot; hold on;
180 173
             plot(  dotC(:,1),dotC(:,2),'w+','MarkerSize',14 );
@@ -219,7 +212,7 @@ disp('loading data...');
219 212
         mRNA_den = zeros(1,Nnucs);  % store densities of mRNA per cell
220 213
         nuc_area = zeros(1,Nnucs); 
221 214
         if showim == 1
222  
-            Plot_mRNA = single(NucLabeled);
  215
+            Plot_mRNA = single(NucLabel);
223 216
         end
224 217
         for i=1:Nnucs; % i = 4
225 218
             nn = Nucs_list(i);
8  jacquestiffread.m
@@ -131,6 +131,8 @@
131 131
         switch TIF.entry_tag
132 132
             case 254
133 133
                 TIF.NewSubfiletype = entry.val;
  134
+                
  135
+                           
134 136
             case 256         % image width - number of column
135 137
                 IMG.width          = entry.val;
136 138
                 
@@ -212,7 +214,7 @@
212 214
 
213 215
 Data.LSM_info=LSM_info;
214 216
 
215  
-numberofstacks = round((count-2)/Data.LSM_info.DimensionZ);
  217
+numberofstacks = (count-1)/Data.LSM_info.DimensionZ;
216 218
 
217 219
 for i =1:numberofstacks
218 220
     for j=1:Data.LSM_info.DimensionZ
@@ -231,7 +233,7 @@
231 233
 save(filenameout,'Datas');
232 234
 disp('data saved'); 
233 235
 
234  
-fclose(TIF.file)
  236
+fclose(TIF.file);
235 237
 
236 238
 end
237 239
 
@@ -336,7 +338,7 @@
336 338
            
337 339
            end
338 340
            
339  
-           entry.val(i)= entry.val(i)+TEMPDat.counterr*(2^32-1);
  341
+           entry.val(i)= entry.val(i)+TEMPDat.counterr*(2^32);
340 342
            
341 343
            TEMPDat.pastposi=entryvaltemp(i);
342 344
            
10  jacquestiffread.m~
@@ -131,6 +131,8 @@ while  TIF.img_pos ~= 0
131 131
         switch TIF.entry_tag
132 132
             case 254
133 133
                 TIF.NewSubfiletype = entry.val;
  134
+                
  135
+                           
134 136
             case 256         % image width - number of column
135 137
                 IMG.width          = entry.val;
136 138
                 
@@ -212,7 +214,7 @@ end
212 214
 
213 215
 Data.LSM_info=LSM_info;
214 216
 
215  
-numberofstacks = round((count-2)/Data.LSM_info.DimensionZ);
  217
+numberofstacks = (count-1)/Data.LSM_info.DimensionZ;
216 218
 
217 219
 for i =1:numberofstacks
218 220
     for j=1:Data.LSM_info.DimensionZ
@@ -231,7 +233,7 @@ Datas.TEMPDat = TIF.TEMPDat;
231 233
 save(filenameout,'Datas');
232 234
 disp('data saved'); 
233 235
 
234  
-fclose(TIF.file)
  236
+fclose(TIF.file);
235 237
 
236 238
 end
237 239
 
@@ -342,9 +344,9 @@ else
342 344
            
343 345
            end
344 346
            
  347
+           TEMPDat.entryval = entry.val;   % Added export for troubleshooting
345 348
            TIF.TEMPDat=TEMPDat;
346  
-           Datas.TEMPDat=TEMPDat;
347  
-           
  349
+         
348 350
            
349 351
         
350 352
   else

0 notes on commit 8400ff2

Please sign in to comment.
Something went wrong with that request. Please try again.