diff --git a/matlab/+Batch/batchFeatures.m b/matlab/+Batch/batchFeatures.m
new file mode 100644
index 0000000000000000000000000000000000000000..623c22e5bbfbed6b583edad4f150b1b0b93f3f4d
--- /dev/null
+++ b/matlab/+Batch/batchFeatures.m
@@ -0,0 +1,47 @@
+function batchFeatures(conn,CONSTANTS,nProcessors)
+
+algorithms=Read.getAlgorithms(conn,'feature');
+if isempty(algorithms)
+ return
+end
+segParams=Read.getSegmentationParams(conn);
+CONSTANTS=CellFeatures.resetFeatures(conn,CONSTANTS,algorithms);
+
+if ~exist('nProcessors','var') || isempty(nProcessors)
+ nProcessors=feature('numcores');
+end
+p=gcp('nocreate');
+if isempty(p) || p.NumWorkers~=nProcessors
+ delete(p);
+ p=[];
+end
+if isempty(p)
+ p=parpool(nProcessors);
+end
+
+strDB=conn.DataSource;
+for iAlgorithm=1:length(algorithms)
+ tNext=1;
+ for idp=1:ceil(CONSTANTS.imageData.NumberOfFrames/p.NumWorkers)
+ featureTable=Composite();
+ spmd
+ tTarget=tNext+labindex-1; % labIndex is 1 based
+ if tTarget<=CONSTANTS.imageData.NumberOfFrames
+ cx=openDB(strDB);
+ featureTable=algorithms(iAlgorithm).function(cx,tTarget,CONSTANTS,algorithms(iAlgorithm).params,segParams);
+ close(cx);
+ else
+ featureTable=[];
+ end
+ end
+ for i=1:length(featureTable)
+ if ~isempty(featureTable{i})
+ CellFeatures.writeFeatureTable(conn,featureTable{i});
+ end
+ end
+ tNext=tNext+p.NumWorkers;
+ end
+
+end
+
+4;
\ No newline at end of file
diff --git a/matlab/+Batch/processMovie.m b/matlab/+Batch/processMovie.m
index 5e342bef44bb4959cb98abb0085b8c7495c9d8e1..a2ccae63f9d6677cce4f03a1135b15608f584021 100644
--- a/matlab/+Batch/processMovie.m
+++ b/matlab/+Batch/processMovie.m
@@ -22,5 +22,8 @@ tx=tic();
Smooth.Classify.goClassifyMitosis(conn,CONSTANTS);
Smooth.extFamily(conn);
tProcess.tClassify=toc(tx);
+tx=tic();
+Batch.batchFeatures(conn,CONSTANTS);
+tProcess.tFeatures=toc();
close(conn)
diff --git a/matlab/+CellFeatures/cytoToNuclearRatio.m b/matlab/+CellFeatures/cytoToNuclearRatio.m
index c589c19762cbf643481ef4fc3660ab1bb159f137..a80c2d40d3489f43eab67efdb9f18b9b13571b4e 100644
--- a/matlab/+CellFeatures/cytoToNuclearRatio.m
+++ b/matlab/+CellFeatures/cytoToNuclearRatio.m
@@ -1,105 +1,57 @@
-function cytoToNuclearRatio(conn,CONSTANTS,args,t)
+function tblCellFeatures=cytoToNuclearRatio(conn,t,CONSTANTS,featureParams,segParams)
+
+tblCellFeatures=table();
params.channel=1;
params.name='CytoToNuclearRatio';
params.color="#00ffff";
% check for param override
-params=Segment.getParams(params,args);
-
-% segmentation params -- for min radius
-algs=Read.getAlgorithms(conn,'segment');
-args=algs.params;
-segParams=Segment.getDefaultSegParams();
-segParams=Segment.getParams(segParams,args);
-
-rgRadius=Ensemble.getEnsembleRadius(segParams.minimumRadius_um);
-resolution_um=CONSTANTS.imageData.PixelPhysicalSize; % um per pixel
-min_radius_pixels = rgRadius ./ resolution_um(1);
-
-% todo -- handle anistropy in 3d!
+featureParams=Segment.getParams(params,featureParams);
-cmd=['select cellID from tblCells where time=' num2str(t)];
-q=ljsFetch(conn,cmd);
-cellIDs=cell2mat(q);
-if isempty(cellIDs)
- return
+cellFeatures=[];
+cellList=Read.getCellsTime(conn,t);
+if isempty(cellList)
+ return;
end
+cellList=Read.setCellsIdxPts(cellList,CONSTANTS);
+
im = MicroscopeData.Reader('imageData',CONSTANTS.imageData, 'chanList',...
- params.channel,'timeRange',[t t], 'outType','single','prompt',false);
+ featureParams.channel,'timeRange',[t t], 'outType','single','prompt',false);
if all(isnan(im(:)))
return;
end
is3d=size(im,3)>1;
labelCells=0*im;
dInterior=labelCells;
-for i=1:length(cellIDs)
- if is3d
- c=Read.getCell3D(conn,cellIDs(i));
- idx=sub2ind(size(im),c.pts(:,2),c.pts(:,1),c.pts(:,3));
- else
- c=Read.getCell(conn,cellIDs(i));
- idx=sub2ind(size(im),c.pts(:,2),c.pts(:,1));
- end
+for i=1:length(cellList)
+ idx=cellList(i).idxPts;
bwInterior=logical(0*im);
bwInterior(idx)=1;
+ % we need bwdist for just this cell
bwd=bwdist(~bwInterior);
+ % isolated bwdist put in combo image, indexed by labelCells
dInterior=dInterior+bwd;
- labelCells(idx)=c.cellID;
+ labelCells(idx)=cellList(i).cellID;
end
% find region surrounding segmentation
[dExterior,idxd]=bwdist(labelCells);
-
-for i=1:length(rgRadius)
- li=labelCells;
- lo=labelCells(idxd);
- radius=rgRadius(i);
- lo(dExterior>radius)=0;
- % should this be 0, or min_radius_pixels? i.e. buffer between?
- lo(dInterior>=radius)=0;
- li(dInterior<radius)=0;
-
- labelOuter(:,:,:,i)=lo;
- labelInner(:,:,:,i)=li;
-end
-
-for i=1:length(cellIDs)
- if is3d
- c=Read.getCell3D(conn,cellIDs(i));
- else
- c=Read.getCell(conn,cellIDs(i));
- end
- if isfield(c,'segCC')
- iRadius=c.segCC-floor(c.segCC);
- if 0==iRadius
- iRadius=1;
- else
- prec=ceil(log10(length(rgRadius)+1));% # of digits in the radius tag
- iRadius=iRadius*10^prec;
- iRadius=int32(iRadius);
- end
- else
- iRadius=1;
- end
+%
+for i=1:length(cellList)
- idxInner=find(labelInner(:,:,:,iRadius)==c.cellID);
- idxOuter=find(labelOuter(:,:,:,iRadius)==c.cellID);
+ % interior pixels > 2 pixels from cell surface
+ idxInner = find( labelCells==cellList(i).cellID & dInterior>=2);
+ % pixels <=2 pixels from cell surface
+ idxInteriorBoundary=find( labelCells==cellList(i).cellID & dInterior<2);
+ % exterior pixels nearest to this cell
+ idxExterior=idxd(find( labelCells(idxd)==cellList(i).cellID & dExterior<=3 & dExterior>0));
+ idxOuter=[idxInteriorBoundary;idxExterior];
ratio=getErkRatio(im,idxOuter,idxInner);
- if isnan(ratio)
- cmd=['delete from tblCellFeatures where cellID=' num2str(cellIDs(i))];
- exec(conn,cmd);
- else
- ratio=round(ratio,4);
- cmd=['insert or ignore into tblCellFeatures(cellID) values (' ...
- num2str(cellIDs(i)) ')'];
- exec(conn,cmd);
- cmd=['update tblCellFeatures set ' params.name '=' num2str(ratio) ...
- ' where cellID=' num2str(cellIDs(i))];
- exec(conn,cmd);
- end
-
+ cellFeatures=[cellFeatures;cellList(i).cellID,ratio];
end
-
+tblCellFeatures.cellID=cellFeatures(:,1);
+tblCellFeatures.(featureParams.name)=cellFeatures(:,2);
+4;
function ratio=getErkRatio(im,idxOuter,idxInner)
ratio=NaN;
diff --git a/matlab/+CellFeatures/netIntensity.m b/matlab/+CellFeatures/netIntensity.m
index a8cff5b68f015a486288f4fe04fb8b1bd55b4d21..528fc86a1c67576c1a881e3c02e49bf871d8871b 100644
--- a/matlab/+CellFeatures/netIntensity.m
+++ b/matlab/+CellFeatures/netIntensity.m
@@ -1,30 +1,31 @@
-function netIntensity(conn,CONSTANTS,args,t)
+
+function tblCellFeatures=netIntensity(conn,t,CONSTANTS,featureParams,segParams)
+
+tblCellFeatures=table();
+cellFeatures=[];
params.channel=1;
params.name='NetIntensity';
params.color="#ff0000";
% check for param override
-params=Segment.getParams(params,args);
+featureParams=Segment.getParams(params,featureParams);
-cmd=['select cellID from tblCells where time=' num2str(t)];
-q=ljsFetch(conn,cmd);
-cellIDs=cell2mat(q);
+cellList=Read.getCellsTime(conn,t);
+if isempty(cellList)
+ return;
+end
+cellList=Read.setCellsIdxPts(cellList,CONSTANTS);
im = MicroscopeData.Reader('imageData',CONSTANTS.imageData, 'chanList',...
- params.channel,'timeRange',[t t], 'outType','single','prompt',false);
+ featureParams.channel,'timeRange',[t t], 'outType','single','prompt',false);
if all(im(:)==0)
return
end
bias=median(im(:));
-for i=1:length(cellIDs)
- c=Read.getCell(conn,cellIDs(i));
- idxSegmentation=sub2ind(size(im),c.pts(:,2),c.pts(:,1));
-
- net=sum(im(idxSegmentation)-bias);
- cmd=['insert or ignore into tblCellFeatures(cellID) values (' ...
- num2str(cellIDs(i)) ')'];
- exec(conn,cmd);
- cmd=['update tblCellFeatures set ' params.name '=' num2str(net) ...
- ' where cellID=' num2str(cellIDs(i))];
- exec(conn,cmd);
+for i=1:length(cellList)
+ net=sum(im(cellList(i).idxPts)-bias);
+ cellFeatures=[cellFeatures;cellList(i).cellID,net];
end
+tblCellFeatures.cellID=cellFeatures(:,1);
+tblCellFeatures.(featureParams.name)=cellFeatures(:,2);
+4;
diff --git a/matlab/+CellFeatures/netNuclearEnvelope.m b/matlab/+CellFeatures/netNuclearEnvelope.m
index 1f249b1e04b86df658fd44ad833e7631daf2a1bd..fb26140c95958aa39cd7c00b8dab8c75fe869033 100644
--- a/matlab/+CellFeatures/netNuclearEnvelope.m
+++ b/matlab/+CellFeatures/netNuclearEnvelope.m
@@ -1,109 +1,49 @@
-function netNuclearEnvelope(conn,CONSTANTS,args,t)
+
+function tblCellFeatures=netNuclearEnvelope(conn,t,CONSTANTS,featureParams,segParams)
+
+tblCellFeatures=table();
params.channel=1;
params.name='netNuclearEnvelope';
params.color="#ff0000";
% check for param override
-params=Segment.getParams(params,args);
-
-% segmentation params -- for min radius
-algs=Read.getAlgorithms(conn,'segment');
-args=algs.params;
-segParams=Segment.getDefaultSegParams();
-segParams=Segment.getParams(segParams,args);
-
-rgRadius=Ensemble.getEnsembleRadius(segParams.minimumRadius_um);
-resolution_um=CONSTANTS.imageData.PixelPhysicalSize; % um per pixel
-min_radius_pixels = rgRadius ./ resolution_um(1);
+featureParams=Segment.getParams(params,featureParams);
-% todo -- handle anistropy in 3d!
-
-cmd=['select cellID from tblCells where time=' num2str(t)];
-q=ljsFetch(conn,cmd);
-cellIDs=cell2mat(q);
+cellList=Read.getCellsTime(conn,t);
+if isempty(cellList)
+ return;
+end
+cellList=Read.setCellsIdxPts(cellList,CONSTANTS);
im = MicroscopeData.Reader('imageData',CONSTANTS.imageData, 'chanList',...
- params.channel,'timeRange',[t t], 'outType','single','prompt',false);
+ featureParams.channel,'timeRange',[t t], 'outType','single','prompt',false);
if any(isnan(im))
return;
end
-is3d=size(im,3)>1;
labelCells=0*im;
dInterior=labelCells;
-for i=1:length(cellIDs)
- if is3d
- c=Read.getCell3D(conn,cellIDs(i));
- idx=sub2ind(size(im),c.pts(:,2),c.pts(:,1),c.pts(:,3));
- else
- c=Read.getCell(conn,cellIDs(i));
- idx=sub2ind(size(im),c.pts(:,2),c.pts(:,1));
- end
+% get separate distance transform for each cell
+for i=1:length(cellIDs)
bwInterior=logical(0*im);
- bwInterior(idx)=1;
+ bwInterior(cellList(i).idxPts)=1;
bwd=bwdist(~bwInterior);
dInterior=dInterior+bwd;
- labelCells(idx)=c.cellID;
+ labelCells(idx)=cellList(i).cellID;
end
% find region surrounding segmentation
[dExterior,idxd]=bwdist(labelCells);
-for i=1:length(rgRadius)
- li=labelCells;
- lo=labelCells(idxd);
- radius=rgRadius(i);
- lo(dExterior>radius)=0;
- % should this be 0, or min_radius_pixels? i.e. buffer between?
- lo(dInterior>=radius)=0;
- li(dInterior>=radius)=0;
-
- labelOuter(:,:,:,i)=lo;
- labelInner(:,:,:,i)=li;
-end
bias=median(im(:));
for i=1:length(cellIDs)
- if is3d
- c=Read.getCell3D(conn,cellIDs(i));
- else
- c=Read.getCell(conn,cellIDs(i));
- end
- if isfield(c,'segCC')
- iRadius=c.segCC-floor(c.segCC);
- if 0==iRadius
- iRadius=1;
- else
- prec=ceil(log10(length(rgRadius)+1));% # of digits in the radius tag
- iRadius=iRadius*10^prec;
- iRadius=int32(iRadius);
- end
- else
- iRadius=1;
- end
-
- idxInner=find(labelInner(:,:,:,iRadius)==c.cellID);
- idxOuter=find(labelOuter(:,:,:,iRadius)==c.cellID);
- net=mean( [im(idxInner)-bias;im(idxOuter)-bias]);
-
- net=round(net,4);
- cmd=['insert or ignore into tblCellFeatures(cellID) values (' ...
- num2str(cellIDs(i)) ')'];
- exec(conn,cmd);
- cmd=['update tblCellFeatures set ' params.name '=' num2str(net) ...
- ' where cellID=' num2str(cellIDs(i))];
- exec(conn,cmd);
-end
-
-function ratio=getErkRatio(im,idxOuter,idxInner)
-ratio=NaN;
-imOuter=im(idxOuter);
-imOuter(0==imOuter)=[];
-if length(imOuter)<5
- return
-end
-imInner=im(idxInner);
-imInner(0==imInner)=[];
-if length(imInner)/length(idxInner)<0.1
- return
+ % read signal within 2 pixels on either side of the boundary
+ idxInner = find( labelCells==cellList(i).cellID & dInterior<=2);
+ idxOuter = idxd(find( labelCells(idxd)==cellList(i).cellID & dExterior<=2 & dExterior>0));
+ net=mean( [im(idxInner)-bias;im(idxOuter)-bias]);
+ net=round(net,4);
+ cellFeatures=[cellFeatures;cellList(i).cellID,net];
end
-ratio=median(imOuter)/median(imInner);
+tblCellFeatures.cellID=cellFeatures(:,1);
+tblCellFeatures.(featureParams.name)=cellFeatures(:,2);
+4;
diff --git a/matlab/+CellFeatures/writeFeatureTable.m b/matlab/+CellFeatures/writeFeatureTable.m
new file mode 100644
index 0000000000000000000000000000000000000000..5183cd0ac96a4e932f054f482a790254e859f01f
--- /dev/null
+++ b/matlab/+CellFeatures/writeFeatureTable.m
@@ -0,0 +1,42 @@
+% should be cellID and 1 feature. write to db.
+function writeFeatureTable(conn,featureTable)
+
+if size(featureTable,2)~=2 || ~strcmp(featureTable.Properties.VariableNames{1},'cellID')
+ fprintf(2,'writeFeatureTable - invalid table format\n');
+ return;
+end
+
+% if isnan(ratio)
+% cmd=['delete from tblCellFeatures where cellID=' num2str(cellIDs(i))];
+% exec(conn,cmd);
+% else
+% ratio=round(ratio,4);
+% cmd=['insert or ignore into tblCellFeatures(cellID) values (' ...
+% num2str(cellIDs(i)) ')'];
+% exec(conn,cmd);
+% cmd=['update tblCellFeatures set ' params.name '=' num2str(ratio) ...
+% ' where cellID=' num2str(cellIDs(i))];
+% exec(conn,cmd);
+% end
+cmdDelete='delete from tblCellFeatures where cellID=?';
+psDelete=conn.handle.prepareStatement(cmdDelete);
+cmdInsert='insert or ignore into tblCellFeatures(cellID) values (?)';
+psInsert=conn.handle.prepareStatement(cmdInsert);
+cmdUpdate=['update tblCellFeatures set ' featureTable.Properties.VariableNames{2}...
+ '=? where cellID=?'];
+psUpdate=conn.handle.prepareStatement(cmdUpdate);
+
+for i=1:size(featureTable,1)
+ if isnan(featureTable{i,2})
+ psDelete.setDouble(1,featureTable{i,1});
+ Write.tryUpdateDB(psDelete);
+ else
+ psInsert.setDouble(1,featureTable{i,1});
+ Write.tryUpdateDB(psInsert);
+ % set feature value
+ psUpdate.setDouble(1,featureTable{i,2});
+ % set cellID
+ psUpdate.setDouble(2,featureTable{i,1});
+ Write.tryUpdateDB(psUpdate);
+ end
+end
diff --git a/matlab/+Segment/thresholdImages.m b/matlab/+Segment/thresholdImages.m
index f094a6a05c4aa2bd8c1165326ddb878f0c5be16a..95c1994879f0b61677bd1c03eb9d62e18eae994b 100644
--- a/matlab/+Segment/thresholdImages.m
+++ b/matlab/+Segment/thresholdImages.m
@@ -107,6 +107,7 @@ else
pixelIdxList={rp.PixelIdxList};
end
idxReset=find(convDef<1.1);
-pixelIdxReset=vertcat(pixelIdxList{idxReset});
-bwLoG(pixelIdxReset)=0;
-
+if ~isempty(idxReset)
+ pixelIdxReset=vertcat(pixelIdxList{idxReset});
+ bwLoG(pixelIdxReset)=0;
+end
diff --git a/matlab/+Smooth/+Classify/goClassifyMitosis.m b/matlab/+Smooth/+Classify/goClassifyMitosis.m
index a1180366c2c1424dd2d23c7cb0d9b7c43a134b0f..81e8f5d7b5aba99584f172bb9ed976b51a463b50 100644
--- a/matlab/+Smooth/+Classify/goClassifyMitosis.m
+++ b/matlab/+Smooth/+Classify/goClassifyMitosis.m
@@ -2,6 +2,11 @@
% pass in [conn,CONSTANTS] (from openDB) or strDB (full path to .LEVER)
function classifyMitosis(conn,CONSTANTS)
+reseg=Read.getAlgorithms(conn,'reseg');
+if isempty(reseg) || ~strcmp(reseg.params.mitosisDetection,'true')
+ return;
+end
+
mParams=Read.getMitosisParams(conn);
targetTracks=Smooth.getTargetTracks(conn,5);
strDB=conn.DataSource;