From 816c4127aab108e19e2136fcdf1f894ffbfe95dc Mon Sep 17 00:00:00 2001 From: Anna Yoney Date: Tue, 15 Apr 2014 11:48:53 -0400 Subject: [PATCH 1/3] Test commit --- Cytoo/colonyColorPointPlot.m | 17 ++- Cytoo/getDimsFromLogFile.m | 82 +++++++++++ Cytoo/mkFullCytooPlotPeaks.m | 44 ++++-- Cytoo/peaksToColonies.m | 7 +- Cytoo/runTileLoop.m | 2 +- paramFiles/setUserParam20XAnna.m | 226 +++++++++++++++++++++++++++++++ runFiles/runDeltaVisionTile.m | 136 +++++++++++++++++++ runFiles/runFullTile.m | 3 +- 8 files changed, 494 insertions(+), 23 deletions(-) create mode 100644 Cytoo/getDimsFromLogFile.m create mode 100755 paramFiles/setUserParam20XAnna.m create mode 100644 runFiles/runDeltaVisionTile.m diff --git a/Cytoo/colonyColorPointPlot.m b/Cytoo/colonyColorPointPlot.m index 5897a72..0ac31b8 100644 --- a/Cytoo/colonyColorPointPlot.m +++ b/Cytoo/colonyColorPointPlot.m @@ -1,4 +1,4 @@ -function colonyColorPointPlot(col,dcols,ps,climits,newfig,rescale_fac,centerpos) +function colonyColorPointPlot(col,dcols,ps,climits,newfig,rescale_fac,centerpos,dofilter) %function colonyColorPointPlot(col,dcols) %--------------------------------------- %Colony scatter plot with points colored by data @@ -7,6 +7,7 @@ function colonyColorPointPlot(col,dcols,ps,climits,newfig,rescale_fac,centerpos) % if length==2, use 2nd one for norm %ps = pointsize (default 12); + xdat=col.data(:,1); ydat=col.data(:,2); @@ -37,8 +38,12 @@ function colonyColorPointPlot(col,dcols,ps,climits,newfig,rescale_fac,centerpos) centerpos=1; end +if ~exist('dofilter','var') + dofilter=1; +end + if newfig -figure; + figure; end if centerpos @@ -46,12 +51,16 @@ function colonyColorPointPlot(col,dcols,ps,climits,newfig,rescale_fac,centerpos) ydat=bsxfun(@minus,ydat,mean(ydat)); end -includeinds=sqrt(xdat.*xdat+ydat.*ydat) < 30*50+5; +if dofilter + includeinds=sqrt(xdat.*xdat+ydat.*ydat) < 30*50+5; +else + includeinds=1:length(xdat); +end xdat=xdat(includeinds); ydat=ydat(includeinds); coldat=coldat(includeinds); colormap('jet'); -scatter(rescale_fac*xdat,rescale_fac*ydat,ps,coldat); +scatter(rescale_fac*xdat,rescale_fac*ydat,ps,coldat,'filled'); minx=min(xdat); maxx=max(xdat); miny=min(ydat); maxy=max(ydat); diff --git a/Cytoo/getDimsFromLogFile.m b/Cytoo/getDimsFromLogFile.m new file mode 100644 index 0000000..cd9b2c6 --- /dev/null +++ b/Cytoo/getDimsFromLogFile.m @@ -0,0 +1,82 @@ +function [dims, wavenames]=getDimsFromLogFile(direc) +%[dims wavenames]=getDimsFromLogFile(direc) +%----------------------------------------------- +%Function to find dimensions of the tiling and +%the names of the wavelengths from the .dv.log file outputted by +%deltavision +%direc = output directory containing the .dv.log file +% assumes only one .dv.log file in that directory. +%dims = 2 component vector containing dimension +%wavenames = cell array containing wavelength names + +verbose = 0; + +fnames=dir([direc filesep '*log']); +scanfile=[direc filesep fnames(1).name]; +ff=fopen(scanfile); + +strtofind1='DO '; +strtofind2='CHANNEL'; +strtofind3='Stage coordinates:'; +strtofind4='Image '; + +tline=fgetl(ff); +q=1; +Y=1; + +while ischar(tline) + + k1=strfind(tline,strtofind1); + k2=strfind(tline,strtofind2); + k4=strfind(tline,strtofind3); + n=strfind(tline,strtofind4); + + if k1 + if verbose + disp(tline); + end + panels=str2double(tline(k1+3:end)); + end + + if k2 + if verbose + disp(tline); + end + + k3=strfind(tline,','); + wavenames{q}=tline(k3(1)+1:k3(2)-1); + q=q+1; + end + + if n + n1=strfind(tline,'.'); + imageNum=str2double(tline(n+6:n1-1)); + end + + if k4 + if verbose + disp(tline); + end + + k5=strfind(tline,','); + y=str2double(tline(k5(1)+1:k5(2)-1)); + + if Y==1 + y1=y; + Y=2; + end + + if Y==2 && y~=y1 + cols = (imageNum-1)./length(wavenames); + Y=3; + end + + end + + tline=fgetl(ff); + +end + +rows = floor(panels./cols); + +dims = [rows,cols]; diff --git a/Cytoo/mkFullCytooPlotPeaks.m b/Cytoo/mkFullCytooPlotPeaks.m index 49d49ff..7b2bc5d 100644 --- a/Cytoo/mkFullCytooPlotPeaks.m +++ b/Cytoo/mkFullCytooPlotPeaks.m @@ -13,25 +13,41 @@ pp=load(matfile); peaks=pp.peaks; ac=pp.acoords; +col=pp.plate1.colonies; + +totcells = 0; +for ii=1:length(peaks) + totcells=totcells+size(peaks{ii},1); +end + +ncolumn=size(col(1).data,2); +alldat=zeros(totcells,ncolumn); + + cc=colorcube(19); figure; hold on; +q=1; + for ii=1:length(peaks) if ~isempty(peaks{ii}) - toadd=[ac(ii).absinds(2) ac(ii).absinds(1)]; - %toadd=[0 0]; - dtoplot=bsxfun(@plus,peaks{ii}(:,1:2),toadd); - if returndat - alldat(q:(q+col(ii).ncells-1),:)=[dtoplot col(ii).data(:,3:end)]; - q=q+col(ii).ncells; - end - plot(dtoplot(:,2),dtoplot(:,1),'.','Color',cc(mod(ii,19)+1,:)); -% txtcolor=cc(mod(ii+10,19)+1,:); -% cen=[col(ii).center(2)+toadd(2) col(ii).center(1)+toadd(1)]; -% if printnum && col(ii).ncells > 200 -% text(cen(1),cen(2),int2str(ii),'Color','m'); -% end - %disp(int2str(ii)); + toadd=[ac(ii).absinds(2) ac(ii).absinds(1)]; + %toadd=[0 0]; + dtoplot=bsxfun(@plus,peaks{ii}(:,1:2),toadd); + if returndat + alldat(q:(q+size(peaks{ii},1)-1),:)=[dtoplot peaks{ii}(:,3:end)]; + q=q+size(peaks{ii},1); + + end + + text(mean(dtoplot(:,2)),mean(dtoplot(:,1)),int2str(ii),'Color','k'); + plot(dtoplot(:,2),dtoplot(:,1),'.','Color',cc(mod(ii,19)+1,:)); + % txtcolor=cc(mod(ii+10,19)+1,:); + % cen=[col(ii).center(2)+toadd(2) col(ii).center(1)+toadd(1)]; + % if printnum && col(ii).ncells > 200 + % text(cen(1),cen(2),int2str(ii),'Color','m'); + % end + %disp(int2str(ii)); end end diff --git a/Cytoo/peaksToColonies.m b/Cytoo/peaksToColonies.m index 4320c35..1f42996 100644 --- a/Cytoo/peaksToColonies.m +++ b/Cytoo/peaksToColonies.m @@ -1,9 +1,10 @@ -function [colonies peaks]=peaksToColonies(matfile) +function [colonies, peaks]=peaksToColonies(matfile) -pp=load(matfile,'peaks','acoords','imgfiles','dims'); +pp=load(matfile,'peaks','acoords','imgfiles','dims','si'); peaks=pp.peaks; ac=pp.acoords; dims=pp.dims; +si=pp.si; peaks=removeDuplicateCells(peaks,ac); @@ -40,7 +41,7 @@ %Make colony structure for ii=1:length(groups) cellstouse=allinds==ii; - colonies(ii)=colony(alldat(cellstouse,:),ac,dims,[],pp.imgfiles); + colonies(ii)=colony(alldat(cellstouse,:),ac,dims,si,pp.imgfiles); end %put data back into peaks diff --git a/Cytoo/runTileLoop.m b/Cytoo/runTileLoop.m index c6daf13..406c8bb 100644 --- a/Cytoo/runTileLoop.m +++ b/Cytoo/runTileLoop.m @@ -4,7 +4,7 @@ function runTileLoop(direc,chans,imgsperprocessor,nloop,maxims,bIms,nIms,paramfi if s > 0 matlabpool close; end -matlabpool('local',12); +matlabpool('local'); parfor ii=1:nloop n1=(ii-1)*imgsperprocessor+1; n2=min(ii*imgsperprocessor,maxims); diff --git a/paramFiles/setUserParam20XAnna.m b/paramFiles/setUserParam20XAnna.m new file mode 100755 index 0000000..88538ef --- /dev/null +++ b/paramFiles/setUserParam20XAnna.m @@ -0,0 +1,226 @@ +function setUserParam20XAnna(img) +% +% setUserParamsCCC20x(img) +% +% Contains master set of comments on how to adjust parameters and other +% hints. Parameters based on 110118 images. The argument img is used only +% to extract dimensions, use img=[] to use default +% +% Note on file compression: to go from 1024x1344 1.6mb .png to ~110kb +% 512x672, mask both nuc and cyto images with dilation of the total cell +% mask (ie set to ==0 pixels outside of the mask) and +% then save as imwrite(... 'name.jpeg', 'Bitdepth', 16) and use default +% compression. (might try setting backgnd to min(img) and using 8 bit) + +global userParam + +rescale = 1.0; +if nargin && ~isempty(img) && (min(size(img)) > 512) + % do not use, + %rescale = 0.5; +end + +fprintf(1, 'setUserParam10XmSCEDS(): called to define params, rescale factor= %d\n', rescale); + +% When verbose=1 set, image of field of cells produced with diagnostics. If +% newFigure=1 these will pile up for successive times and eventually crash MATLAB +% because of memory limitations. Either run in debug mode and kill by hand or set +% newFigure=0. +userParam.newFigure = 0; +userParam.errorStr = []; + +%%%%%%%%%%%%%%% used in segmentCells() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +userParam.verboseSegmentCells = 0; + +% use at least one of these methods to find nuclei. +% findNucThresh() thresholds on a global intensity value determined via +% gauss_thresh(), and then tries to segment composite nuclei. Can miss +% faint nuclei, but these picked up by subsequent call to countNucCtr. If +% nuclei overly large, increase params for gaussian_thresh() below. +% countNucCtr() begins by finding all local max intensity and then +% expands each by looking at contour height of max intensity gradient. Can +% fail to separate nuclei without clear saddle in intensity inbetween, but +% the size of nuclei that are segmented looks very good. +userParam.findNucThresh = 0; +userParam.countNucCtr = 1; + +% filter all 1024x1300 sized images with this gaussian filter to give clean +% local max, and edge detection, 2 will not work +userParam.gaussFilterRadius = 10; + +% for more control over the limits of nuclei, can use a different filter for +% the image that gets passed to edge_thresh_nuc() +% userParam.gaussFilterRadiusEdge = 1; + +% for images with few nuclei, limit size of voronoi polyhedra to max nuc +% radius around center (inferred from max area parameter) +userParam.limitVoronoi = 0; + +% parameter for edge detection. Use 'canny' method in edge() unless get +% nuclei overly large. Check method by on gaussian filtered image +% edge(red, 'canny') vs edge(red). If loosing nuclei use Canny. Set to zero +% to run default. +userParam.useCanny = 1; + +% define threshold for being in cell by two criterion: +% percent of nuclear area in cells > percNucInCell AND +% area of cells > cyto2NucArea * total_nuc_area +% If get huge area defined as 'cell' decrease percNucInCell to miss a few nuclei +% and get better delineation of cell. Suspect large scale background +% variation in green channel. +userParam.percNucInCell = 0.997; +%%%% following used in routine, cell_thresh() that is commented out. +% userParam.cyto2NucArea = 5; + +% If use edgeThreshCyto() to define cytoplasm for each nuclei separately, +% do not need previous 2 parameters. Following verbose plots cyto - backgnd +userParam.verboseEdgeThreshCyto = 0; + +% backgndMethod = 0 assume just the lowest few percent of pixels background (eg +% use with frog AC) +% backgndMethod = 1 define globally from histogram for entire image +% backgndMethod = 2 define locally for each cell based on following.. +% Algorithm looks for gradient of smoothed green, chose a threshold of median +% intensity of points where a gradient is detected. +% or if no gradients AND the max green in V-poly is +% > bckgnd + sclCytoStd*stdb (ie background + some mulitple of std of backgnd) +% then define threshold = cytoHalfMax*(max_img_Vpoly - bckgnd) + bckgnd +% The background can either be defined locally for each Vpoly as imopen(img, big-box)) or as +% one number for entire image. Std can be computed locally but always >= std for entire image +userParam.backgndMethod = -1; +userParam.sclCytoStd = 1; +userParam.cytoHalfMax = 0.5; +userParam.pctLTBckgnd = 0; % for Method=0, percent_less_than_background NB percent thus ~1 +userParam.backdiskrad = 50; %radius for imopen to run on background image before subtracting + +% nucAreaHi used here also but defined below with similar params. + +%%%%%%%%%%%%% Parameters for countNuc(): %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% Used imreconstruct method to find all local max and then filter. +% regmx = imextendedmax(im2, 5) worked somewhat, but gave more false pos. + +userParam.verboseCountNuc = 0; % to print statistics and an image + +userParam.restrictVoronoi=0; % don't allow huge voronoi cells at border + +% Filtering of nuclei done in three steps: +% 1. An optional threshold on allowed min intensity at max. +% 2. A generous test for local max, that should not miss anything and +% also give not too many false +. This is done by imreconstruction +% comparing mask = img + nucIntensityLoc, with dilation of img. Increase +% the nucIntensity local parameter to eliminate false +, decrease if +% missing real nuclei. All local max within a distance of minNucSep are +% merged. +% 3. A filter compares intensity at centroid of each local max with +% intensity in ring defined by radiusMin/Max and demands a contrast of at +% least nucIntensityRange. If loosing real nucs, decrease, if getting false +% positives, increase. +% To adjust the two nucIntensity numbers, run with verbose=1 and look at +% the image. The nuclei after first and second selection shown in different +% colors. If not finding at all obvious nucl, lower thresh in (1) +% +userParam.dontFilterNuc=1; % set to 1 to skip img(center) threshold filtering step +userParam.radiusMin = 5; +userParam.radiusMax = 10; +userParam.minNucSep = 8; +userParam.nucIntensityRange = 2; % value depends on radiusMin/Max +userParam.nucIntensityLoc = 10; % both sould be 5 + +%%%%%%%%%%%%%% Parameters for findNucThresh() | scoreNucShape() %%%%%%%%%%%%%%%%%%%%%%%%%%%% + +userParam.verboseFindNucThresh = 0; + +% locally refine the input mask, produced by global threshold, by making it +% coincide with level defined by max gradient, computed separately for each cc. +userParam.useEdgeThreshNuc = 0; + +% The segmentation method. Watershed is >3x quicker, but needs well defined +% necks in the binary image gotten via thresholding. chenvese will try to +% find the max gradient contour, but can eliminate nuclei of very different +% intensities within same connectec component. +%userParam.segmentation = 'chenvese'; +userParam.segmentation = 'watershed'; +% can refilter sub images destined for segmentation to eliminate internal +% stucture. 0 or [] to skip +userParam.segFilterRadius = 6; + +% try to segment nuclei with > 1 local intensity max if this flag ON +userParam.test1LocalMax = 1; + +%Prior parameters for filtering nuclei based on size/shape, etc from AW +userParam.nucAreaLo =100*rescale^2; +userParam.nucAreaHi = 50000*rescale^2; % not too big +userParam.nucSolidity = 0.50; % get rid of funny shapes +userParam.nucAspectRatio = 3.5; % not too far from circular + +%%%%%%%%%%%%%%% Params for gaussThresh() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Define threshold for grayscale image by fitting gaussian to center of +% intensity distribution and then calling threshold when actual histogram +% counts is > *Excess * model_fit AND intensity > most probable value + +% *Sigma * STD(of model fit, ie ignoring points far in + tail). There is +% buried 'verbose' parameter in this routine. Its assumeed images are +% integer valued, ie not scaled to [0,1] +userParam.gaussThreshExcess = 5; +userParam.gaussThreshSigma = 3; + +%%%%%%%%%%%%%%%% addCellAvr2Stats %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Donut options: compute cytoplasm average from average over a donut shaped region of +% distance RadiusMin to RadiusMax from nucleus that is also inside of Voronoi. +% Units in pixels. +% Set RadiusMax to 0 to skip and integrate over entire cytoplasm. +% RadiusMax > 0 & forceDonut = 0 compute separate stats for donut and cyto +% RadiusMax > 0 & forceDonut = 1 make cyto == donut in all stats. Useful +% when no cyto detected. +% RadiusMax > 0 & intersectDonutCyto = 1, restrict donut to cytoplasm. +% Mutually exclusive with forceDonut=1 +% NB outputData4AWTracker uses only the RadiusMax parameter to decide what data +% to output +% +userParam.donutRadiusMin = 5; % was 5 must be >=0 +userParam.donutRadiusMax = 12; % set to zero to skip +userParam.forceDonut = 1; +userParam.intersectDonutCyto = 0; +if userParam.forceDonut && ~userParam.donutRadiusMax + fprintf(1, 'WARNING inconsistent userParam for donut option.. check\n'); + userParam.forceDonut = 0; +end +if userParam.forceDonut && userParam.intersectDonutCyto + fprintf(1, 'WARNING can not have both forceDonut and intersectDonutCyto ON at same time, forceDonut->0'); + userParam.forceDonut = 0; +end + +% minimum number of pts in cytoplasm (or donut) inorder to count this cell in +% stats (ie give it area, mean etc > 0) +userParam.minPtsCytoplasm = 5; + + + +%%%%%%%%%%%%%%%%%%%%%%%AW added parameters%%%%%%%%%%%%% +% +% +%parameters for reading files. If use sequential is set to 1, then will use +%the strings provided to tracker as a base and append number without zeros +%--i.e. base1.tif, base2.tif etc. If useSequential=0 will use the dir +%command to search for files with the string and ending with suffix +userParam.useSequential=1; +userParam.filesuffix='.TIF'; + + +%Parameters for tracking +userParam.L = 40; % this is the max distance (in pixels) objects in successive + %frames can be separated by and still match +%cost for each process is: +% (cost of process)*(distance to most likely match for process); +% set weights here: +userParam.divisioncost = 1; + +% if cost of all processes exceeds dropcost, +% will not match it. +userParam.dropcost = 50; + +% max number of frames to look ahead: +userParam.maxjump = 3; + +%frame dependent cost (must be a vector of length maxjump) +userParam.skipframecost=1:-0.05:(1-0.05*(userParam.maxjump-1)); diff --git a/runFiles/runDeltaVisionTile.m b/runFiles/runDeltaVisionTile.m new file mode 100644 index 0000000..5981e48 --- /dev/null +++ b/runFiles/runDeltaVisionTile.m @@ -0,0 +1,136 @@ +function runDeltaVisionTile(direc,outfile,paramfile,step,tiffs,prefix) +%runDeltaVisionTile(direc,outfile,paramfile,step,tiffs,prefix) +%--------------------- +%Version of runFullTile for DeltaVision (.dv) files +%Requries MATLAB Bio-Formats files +%Run from directory containing .dv.log file +%For a set of tiled images, runs segmentCells (uses parfor for this), runs +%alignment program for images, outputs in matfile -- peaks -- cell by cells +%list by image, colonies -- colonies data structure +%direc -- image directory +%prefix -- prefix for TIFF images +%outfile -- matfile for output +%step = step to begin at. See code. allows for skipping finding cells etc. +%tiffs = 1 (0) means extract tiff images and save to direc (0 to skip, +%default = 1) + +if ~exist('tiffs','var') + tiffs = 1; +end + +if ~exist('step','var') + step=1; +end + +if ~exist('paramfile','var') + paramfile='setUserParamSC20xIFEDS'; +end + +reader=bfGetReader('path/to/data/file'); +omeMeta=reader.getMetadataStore(); +numChans=omeMeta.getChannelCount(0); +[dims, wavenames]=getDimsFromLogFile('.'); + +for ii = 1:numChans + chanstemp{ii}=['w' char(omeMeta.getChannelEmissionWavelength(0,ii-1))]; +end + +if tiffs == 1 + dvToTiffs(direc,prefix,chanstemp,reader,dims) +end + +chans=orderchans(chanstemp,wavenames); + +ff=folderFilesFromKeyword(direc,chans{1}); +maxims=ff(end-1); + +nloop=12; +imgsperprocessor=ceil(maxims/12); +%generate background image for each channel +if step < 2 + for ii=1:length(chans) + [minI meanI]=mkBackgroundImage(direc,chans{ii},min(500,maxims)); + bIms{ii}=uint16(2^16*minI); + normIm=(meanI-minI); + normIm=normIm.^-1; + normIm=normIm/min(min(normIm)); + nIms{ii}=normIm; + end + + si=size(bIms{1}); + save([direc filesep outfile],'bIms','nIms','dims','si'); +end +%runTileLoop--runs segmentCells in parfor loop, +%send imgsperprocessor to each, nloop = total number necessary +%Assemble Mat Files--puts together matfiles, all data stored as peaks in +%outfile +if step < 3 + load([direc filesep outfile],'bIms','nIms'); + runTileLoop(direc,chans,imgsperprocessor,nloop,maxims,bIms,nIms,paramfile); +end + +%performs a series of pairwise alignments, +%each img is aligned img on top and to the left, pixel overlap +%stored in accords, can also return fully aligned image, but not +%recommended for large numbers of files. +if step < 4 + [acoords]=alignManyPanels(direc,chans{1},1,4,dims,85:150,maxims); + save([direc filesep outfile],'acoords','-append'); +end + +if step < 5 + assembleMatFiles(direc,imgsperprocessor,nloop,outfile); +end +%peaksToColonies generates the colony structure from peaks and accords +%computes alpha volume and then finds all connected components. +if step < 6 + load([direc filesep outfile],'bIms','nIms'); + [colonies, peaks]=peaksToColonies([direc filesep outfile]); + plate1=plate(colonies,dims,direc,chans,bIms,nIms); + save([direc filesep outfile],'plate1','peaks','-append'); +end + +function chans=orderchans(chanstemp,wavenames,nucname) + +if ~exist('nucname','var') + nucname='DAPI'; +end + +kk=strfind(wavenames,nucname); +ii=~cellfun(@isempty,kk); +nuc_ind=find(ii); + +allchans=1:length(wavenames); +nonnucchans=setdiff(allchans,nuc_ind); +chans{1}=chanstemp{nuc_ind}; +for jj=1:length(nonnucchans) + chans{jj+1}=chanstemp{nonnucchans(jj)}; +end + +function dvToTiffs(direc,prefix,chans,reader,dims) + +if ~exist(direc,'file') + mkdir(direc); +end + +for kk=1:length(chans) + for ii=1:dims(1) + for jj=1:dims(2) + if ~mod(ii,2) %even + oldnum=(ii-1)*dims(2)+dims(2)-jj+1; + else + oldnum=(ii-1)*dims(2)+jj; + end + + newnum=(jj-1)*dims(1)+ii; + + reader.setSeries(oldnum-1); + + series_plane = bfGetPlane(reader, kk); + + %t1=imread(fnames(rr==oldnum).name); + newname= [direc filesep prefix '_' chans{kk} '_s' int2str(newnum) '_t1.TIF']; + imwrite(series_plane,newname); + end + end +end diff --git a/runFiles/runFullTile.m b/runFiles/runFullTile.m index 1f0d0ea..2dcde4e 100644 --- a/runFiles/runFullTile.m +++ b/runFiles/runFullTile.m @@ -34,7 +34,8 @@ function runFullTile(direc,outfile,paramfile,step) normIm=normIm/min(min(normIm)); nIms{ii}=normIm; end - save([direc filesep outfile],'bIms','nIms','dims'); + si=size(bIms{1}); + save([direc filesep outfile],'bIms','nIms','dims','si'); end %runTileLoop--runs segmentCells in parfor loop, %send imgsperprocessor to each, nloop = total number necessary From cd1449254dc3a2a5be7aa69660150c40ed263450 Mon Sep 17 00:00:00 2001 From: Anna Yoney Date: Tue, 15 Apr 2014 11:56:34 -0400 Subject: [PATCH 2/3] Updates for DeltaVision files --- Cytoo/runTileLoop.m | 14 -- paramFiles/setUserParam20XAnna.m | 226 ------------------------------- 2 files changed, 240 deletions(-) delete mode 100644 Cytoo/runTileLoop.m delete mode 100755 paramFiles/setUserParam20XAnna.m diff --git a/Cytoo/runTileLoop.m b/Cytoo/runTileLoop.m deleted file mode 100644 index 406c8bb..0000000 --- a/Cytoo/runTileLoop.m +++ /dev/null @@ -1,14 +0,0 @@ -function runTileLoop(direc,chans,imgsperprocessor,nloop,maxims,bIms,nIms,paramfile) - -s=matlabpool('size'); -if s > 0 - matlabpool close; -end -matlabpool('local'); -parfor ii=1:nloop - n1=(ii-1)*imgsperprocessor+1; - n2=min(ii*imgsperprocessor,maxims); - outfile=[direc filesep 'out_' int2str(n1) '.mat']; - runTile(direc,outfile,chans,[n1 n2],bIms,nIms,paramfile); -end -matlabpool close; \ No newline at end of file diff --git a/paramFiles/setUserParam20XAnna.m b/paramFiles/setUserParam20XAnna.m deleted file mode 100755 index 88538ef..0000000 --- a/paramFiles/setUserParam20XAnna.m +++ /dev/null @@ -1,226 +0,0 @@ -function setUserParam20XAnna(img) -% -% setUserParamsCCC20x(img) -% -% Contains master set of comments on how to adjust parameters and other -% hints. Parameters based on 110118 images. The argument img is used only -% to extract dimensions, use img=[] to use default -% -% Note on file compression: to go from 1024x1344 1.6mb .png to ~110kb -% 512x672, mask both nuc and cyto images with dilation of the total cell -% mask (ie set to ==0 pixels outside of the mask) and -% then save as imwrite(... 'name.jpeg', 'Bitdepth', 16) and use default -% compression. (might try setting backgnd to min(img) and using 8 bit) - -global userParam - -rescale = 1.0; -if nargin && ~isempty(img) && (min(size(img)) > 512) - % do not use, - %rescale = 0.5; -end - -fprintf(1, 'setUserParam10XmSCEDS(): called to define params, rescale factor= %d\n', rescale); - -% When verbose=1 set, image of field of cells produced with diagnostics. If -% newFigure=1 these will pile up for successive times and eventually crash MATLAB -% because of memory limitations. Either run in debug mode and kill by hand or set -% newFigure=0. -userParam.newFigure = 0; -userParam.errorStr = []; - -%%%%%%%%%%%%%%% used in segmentCells() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -userParam.verboseSegmentCells = 0; - -% use at least one of these methods to find nuclei. -% findNucThresh() thresholds on a global intensity value determined via -% gauss_thresh(), and then tries to segment composite nuclei. Can miss -% faint nuclei, but these picked up by subsequent call to countNucCtr. If -% nuclei overly large, increase params for gaussian_thresh() below. -% countNucCtr() begins by finding all local max intensity and then -% expands each by looking at contour height of max intensity gradient. Can -% fail to separate nuclei without clear saddle in intensity inbetween, but -% the size of nuclei that are segmented looks very good. -userParam.findNucThresh = 0; -userParam.countNucCtr = 1; - -% filter all 1024x1300 sized images with this gaussian filter to give clean -% local max, and edge detection, 2 will not work -userParam.gaussFilterRadius = 10; - -% for more control over the limits of nuclei, can use a different filter for -% the image that gets passed to edge_thresh_nuc() -% userParam.gaussFilterRadiusEdge = 1; - -% for images with few nuclei, limit size of voronoi polyhedra to max nuc -% radius around center (inferred from max area parameter) -userParam.limitVoronoi = 0; - -% parameter for edge detection. Use 'canny' method in edge() unless get -% nuclei overly large. Check method by on gaussian filtered image -% edge(red, 'canny') vs edge(red). If loosing nuclei use Canny. Set to zero -% to run default. -userParam.useCanny = 1; - -% define threshold for being in cell by two criterion: -% percent of nuclear area in cells > percNucInCell AND -% area of cells > cyto2NucArea * total_nuc_area -% If get huge area defined as 'cell' decrease percNucInCell to miss a few nuclei -% and get better delineation of cell. Suspect large scale background -% variation in green channel. -userParam.percNucInCell = 0.997; -%%%% following used in routine, cell_thresh() that is commented out. -% userParam.cyto2NucArea = 5; - -% If use edgeThreshCyto() to define cytoplasm for each nuclei separately, -% do not need previous 2 parameters. Following verbose plots cyto - backgnd -userParam.verboseEdgeThreshCyto = 0; - -% backgndMethod = 0 assume just the lowest few percent of pixels background (eg -% use with frog AC) -% backgndMethod = 1 define globally from histogram for entire image -% backgndMethod = 2 define locally for each cell based on following.. -% Algorithm looks for gradient of smoothed green, chose a threshold of median -% intensity of points where a gradient is detected. -% or if no gradients AND the max green in V-poly is -% > bckgnd + sclCytoStd*stdb (ie background + some mulitple of std of backgnd) -% then define threshold = cytoHalfMax*(max_img_Vpoly - bckgnd) + bckgnd -% The background can either be defined locally for each Vpoly as imopen(img, big-box)) or as -% one number for entire image. Std can be computed locally but always >= std for entire image -userParam.backgndMethod = -1; -userParam.sclCytoStd = 1; -userParam.cytoHalfMax = 0.5; -userParam.pctLTBckgnd = 0; % for Method=0, percent_less_than_background NB percent thus ~1 -userParam.backdiskrad = 50; %radius for imopen to run on background image before subtracting - -% nucAreaHi used here also but defined below with similar params. - -%%%%%%%%%%%%% Parameters for countNuc(): %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% Used imreconstruct method to find all local max and then filter. -% regmx = imextendedmax(im2, 5) worked somewhat, but gave more false pos. - -userParam.verboseCountNuc = 0; % to print statistics and an image - -userParam.restrictVoronoi=0; % don't allow huge voronoi cells at border - -% Filtering of nuclei done in three steps: -% 1. An optional threshold on allowed min intensity at max. -% 2. A generous test for local max, that should not miss anything and -% also give not too many false +. This is done by imreconstruction -% comparing mask = img + nucIntensityLoc, with dilation of img. Increase -% the nucIntensity local parameter to eliminate false +, decrease if -% missing real nuclei. All local max within a distance of minNucSep are -% merged. -% 3. A filter compares intensity at centroid of each local max with -% intensity in ring defined by radiusMin/Max and demands a contrast of at -% least nucIntensityRange. If loosing real nucs, decrease, if getting false -% positives, increase. -% To adjust the two nucIntensity numbers, run with verbose=1 and look at -% the image. The nuclei after first and second selection shown in different -% colors. If not finding at all obvious nucl, lower thresh in (1) -% -userParam.dontFilterNuc=1; % set to 1 to skip img(center) threshold filtering step -userParam.radiusMin = 5; -userParam.radiusMax = 10; -userParam.minNucSep = 8; -userParam.nucIntensityRange = 2; % value depends on radiusMin/Max -userParam.nucIntensityLoc = 10; % both sould be 5 - -%%%%%%%%%%%%%% Parameters for findNucThresh() | scoreNucShape() %%%%%%%%%%%%%%%%%%%%%%%%%%%% - -userParam.verboseFindNucThresh = 0; - -% locally refine the input mask, produced by global threshold, by making it -% coincide with level defined by max gradient, computed separately for each cc. -userParam.useEdgeThreshNuc = 0; - -% The segmentation method. Watershed is >3x quicker, but needs well defined -% necks in the binary image gotten via thresholding. chenvese will try to -% find the max gradient contour, but can eliminate nuclei of very different -% intensities within same connectec component. -%userParam.segmentation = 'chenvese'; -userParam.segmentation = 'watershed'; -% can refilter sub images destined for segmentation to eliminate internal -% stucture. 0 or [] to skip -userParam.segFilterRadius = 6; - -% try to segment nuclei with > 1 local intensity max if this flag ON -userParam.test1LocalMax = 1; - -%Prior parameters for filtering nuclei based on size/shape, etc from AW -userParam.nucAreaLo =100*rescale^2; -userParam.nucAreaHi = 50000*rescale^2; % not too big -userParam.nucSolidity = 0.50; % get rid of funny shapes -userParam.nucAspectRatio = 3.5; % not too far from circular - -%%%%%%%%%%%%%%% Params for gaussThresh() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Define threshold for grayscale image by fitting gaussian to center of -% intensity distribution and then calling threshold when actual histogram -% counts is > *Excess * model_fit AND intensity > most probable value + -% *Sigma * STD(of model fit, ie ignoring points far in + tail). There is -% buried 'verbose' parameter in this routine. Its assumeed images are -% integer valued, ie not scaled to [0,1] -userParam.gaussThreshExcess = 5; -userParam.gaussThreshSigma = 3; - -%%%%%%%%%%%%%%%% addCellAvr2Stats %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Donut options: compute cytoplasm average from average over a donut shaped region of -% distance RadiusMin to RadiusMax from nucleus that is also inside of Voronoi. -% Units in pixels. -% Set RadiusMax to 0 to skip and integrate over entire cytoplasm. -% RadiusMax > 0 & forceDonut = 0 compute separate stats for donut and cyto -% RadiusMax > 0 & forceDonut = 1 make cyto == donut in all stats. Useful -% when no cyto detected. -% RadiusMax > 0 & intersectDonutCyto = 1, restrict donut to cytoplasm. -% Mutually exclusive with forceDonut=1 -% NB outputData4AWTracker uses only the RadiusMax parameter to decide what data -% to output -% -userParam.donutRadiusMin = 5; % was 5 must be >=0 -userParam.donutRadiusMax = 12; % set to zero to skip -userParam.forceDonut = 1; -userParam.intersectDonutCyto = 0; -if userParam.forceDonut && ~userParam.donutRadiusMax - fprintf(1, 'WARNING inconsistent userParam for donut option.. check\n'); - userParam.forceDonut = 0; -end -if userParam.forceDonut && userParam.intersectDonutCyto - fprintf(1, 'WARNING can not have both forceDonut and intersectDonutCyto ON at same time, forceDonut->0'); - userParam.forceDonut = 0; -end - -% minimum number of pts in cytoplasm (or donut) inorder to count this cell in -% stats (ie give it area, mean etc > 0) -userParam.minPtsCytoplasm = 5; - - - -%%%%%%%%%%%%%%%%%%%%%%%AW added parameters%%%%%%%%%%%%% -% -% -%parameters for reading files. If use sequential is set to 1, then will use -%the strings provided to tracker as a base and append number without zeros -%--i.e. base1.tif, base2.tif etc. If useSequential=0 will use the dir -%command to search for files with the string and ending with suffix -userParam.useSequential=1; -userParam.filesuffix='.TIF'; - - -%Parameters for tracking -userParam.L = 40; % this is the max distance (in pixels) objects in successive - %frames can be separated by and still match -%cost for each process is: -% (cost of process)*(distance to most likely match for process); -% set weights here: -userParam.divisioncost = 1; - -% if cost of all processes exceeds dropcost, -% will not match it. -userParam.dropcost = 50; - -% max number of frames to look ahead: -userParam.maxjump = 3; - -%frame dependent cost (must be a vector of length maxjump) -userParam.skipframecost=1:-0.05:(1-0.05*(userParam.maxjump-1)); From a3146248044662a2a0128e0f1ab589466c8dbbad Mon Sep 17 00:00:00 2001 From: Anna Yoney Date: Tue, 15 Apr 2014 13:28:22 -0400 Subject: [PATCH 3/3] Make compatible with DeltaVision files --- Cytoo/runTileLoop.m | 14 ++++++++++++++ 1 file changed, 14 insertions(+) create mode 100644 Cytoo/runTileLoop.m diff --git a/Cytoo/runTileLoop.m b/Cytoo/runTileLoop.m new file mode 100644 index 0000000..c6daf13 --- /dev/null +++ b/Cytoo/runTileLoop.m @@ -0,0 +1,14 @@ +function runTileLoop(direc,chans,imgsperprocessor,nloop,maxims,bIms,nIms,paramfile) + +s=matlabpool('size'); +if s > 0 + matlabpool close; +end +matlabpool('local',12); +parfor ii=1:nloop + n1=(ii-1)*imgsperprocessor+1; + n2=min(ii*imgsperprocessor,maxims); + outfile=[direc filesep 'out_' int2str(n1) '.mat']; + runTile(direc,outfile,chans,[n1 n2],bIms,nIms,paramfile); +end +matlabpool close; \ No newline at end of file