From 80824eea001ac478921ab215bce721efaac253f4 Mon Sep 17 00:00:00 2001
From: Mikael  BOULLE <mikael.boulle@pasteur.fr>
Date: Mon, 11 Oct 2021 13:57:33 +0200
Subject: [PATCH] Code

---
 Cytokinesis.java | 366 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 366 insertions(+)
 create mode 100644 Cytokinesis.java

diff --git a/Cytokinesis.java b/Cytokinesis.java
new file mode 100644
index 0000000..6036cf3
--- /dev/null
+++ b/Cytokinesis.java
@@ -0,0 +1,366 @@
+package plugins.adufour.hcs;
+
+import java.awt.Color;
+import java.util.Arrays;
+import java.util.List;
+
+import org.apache.poi.ss.usermodel.Workbook;
+
+import icy.math.MathUtil;
+import icy.roi.BooleanMask2D;
+import icy.roi.ROI;
+import icy.roi.ROI2D;
+import icy.roi.ROIUtil;
+import icy.sequence.Sequence;
+import icy.type.collection.CollectionUtil;
+import icy.type.point.Point5D;
+import icy.util.ShapeUtil.BooleanOperator;
+
+import plugins.adufour.activecontours.ActiveContours;
+import plugins.adufour.activecontours.ActiveContours.ROIType;
+import plugins.adufour.blocks.tools.roi.DilateROI;
+import plugins.adufour.ezplug.EzGroup;
+import plugins.adufour.ezplug.EzPlug;
+import plugins.adufour.ezplug.EzStoppable;
+import plugins.adufour.ezplug.EzVarChannel;
+import plugins.adufour.ezplug.EzVarDouble;
+import plugins.adufour.ezplug.EzVarInteger;
+import plugins.adufour.ezplug.EzVarSequence;
+import plugins.adufour.hierarchicalkmeans.HKMeans;
+import plugins.adufour.thresholder.KMeans;
+import plugins.adufour.thresholder.Thresholder;
+import plugins.adufour.workbooks.IcySpreadSheet;
+import plugins.adufour.workbooks.Workbooks;
+import plugins.kernel.roi.descriptor.measure.ROIMassCenterDescriptorsPlugin;
+
+
+public class Cytokinesis extends EzPlug implements EzStoppable
+{
+    
+	
+	
+    @Override
+    public void clean()
+    {
+    }
+    
+    // Define parameters
+    
+    EzVarSequence ezSeq              = new EzVarSequence("Input");
+    
+    EzVarChannel  dapi_channel       = new EzVarChannel("Channel", ezSeq.getVariable(), false);
+    EzVarInteger  dapi_minSize       = new EzVarInteger("Min size", 400, 0, 5000, 50);
+    EzVarInteger  dapi_maxSize       = new EzVarInteger("Max size", 2000, 0, 5000, 50);
+    EzVarDouble   dapi_minIntensity  = new EzVarDouble("Min intensity", 400, 0, 65000, 100);
+   
+    EzVarChannel  tub_channel        = new EzVarChannel("Channel", ezSeq.getVariable(), false);
+    EzVarInteger  tub_classes        = new EzVarInteger("Intensity classes", 2, 2, 100, 1);
+    
+    EzVarChannel  mklp1_channel      = new EzVarChannel("Channel", ezSeq.getVariable(), false);
+    EzVarInteger  mklp1_minSize      = new EzVarInteger("Min size", 30, 0, 1000, 10);
+    EzVarInteger  mklp1_maxSize      = new EzVarInteger("Max size", 100, 0, 1000, 10);
+    EzVarDouble   mklp1_minIntensity = new EzVarDouble("Min intensity", 2000, 0, 65000, 100);
+     
+    EzVarInteger  ring_inflation	 = new EzVarInteger("Ring inflation", 2, 2, 100, 1);
+   
+    
+    
+    @Override
+    protected void initialize()
+    {
+    	
+    	// Create GUI
+    	
+        addEzComponent(ezSeq);
+        
+        addEzComponent(new EzGroup("DAPI", dapi_channel, dapi_minSize, dapi_maxSize, dapi_minIntensity));
+               
+        addEzComponent(new EzGroup("Tubulin", tub_channel, tub_classes));
+        
+        addEzComponent(new EzGroup("MKLP1", mklp1_channel, mklp1_minSize, mklp1_maxSize, mklp1_minIntensity, ring_inflation));            
+    }
+    
+    
+    
+    @Override
+    protected void execute()
+    {
+    	
+    	// Prepare sequence
+    	
+        Sequence seq = ezSeq.getValue(true);
+        
+        seq.removeAllROI();
+        
+        
+        // Prepare collection of statistics  
+        
+        Workbook wb = Workbooks.createEmptyWorkbook();
+        IcySpreadSheet Summary = Workbooks.getSheet(wb, "Summary");
+        IcySpreadSheet FOVs = Workbooks.getSheet(wb, "FOVs");
+        IcySpreadSheet Details = Workbooks.getSheet(wb, "Details");
+        
+        Summary.setRow(0, "Well", "Nb. nuclei", "Nb. midbodies", "Nb. cytokinesis");
+        FOVs.setRow(0, "Well", "Field", "Nb nuclei", "Nb. midbodies", "Nb. cytokinesis");
+        Details.setRow(0, "Well", "Field", "X (dot)", "Y (dot)", "Nb. intersections", "Type");  
+        
+        
+        // Initialise detection loop and total counts
+        
+        int line = 1;
+        int row = 1;
+        
+        int totalNuclei = 0, totalMidbodies = 0, totalCytokinesis = 0;
+        
+	    try  
+	    {
+	        for (int field = 0; field < seq.getSizeT(); field++)
+	        {     
+	        	
+	        	// Prepare FOV
+	        	
+	            Sequence currentField = new Sequence(seq.getImage(field, 0));
+	            
+	            currentField.loadAllData();
+	            currentField.updateChannelsBounds(true);
+	            
+	            
+	            // Follow progression
+	            
+	            getStatus().setCompletion((double) field / seq.getSizeT());
+	            
+	            String prefix = "Frame " + field + ": ";
+	 
+	            if (Thread.currentThread().isInterrupted()) break;
+	            
+	            
+	            // Initialise intermediate counts
+	            
+	            int nucleusNumber = 0;
+	            int midbodyNumber = 0;
+	            int cytokinesisNumber = 0; 
+	            
+	            
+	            // 0. Check for empty images
+	            
+	            if (currentField.getChannelMax(dapi_channel.getValue()) < dapi_minIntensity.getValue())
+	            {
+	            	System.out.println(prefix + "Empty image; intensity in channel " + dapi_channel.getValue() + " < " + dapi_minIntensity.getValue());
+	            	continue;
+	            }  
+	            
+	                 	
+	            // 1. Find DAPI blobs (=> nuclei)
+	                  
+	            	
+	            	// 1.a) HK-Means
+	            
+	            	List<ROI> dapi_roi = HKMeans.hKMeans(currentField, 0, dapi_channel.getValue(), 1.0, (byte) 10,
+	            										 dapi_minSize.getValue(), dapi_maxSize.getValue(),
+	                    			 					 dapi_minIntensity.getValue(), getStatus());
+	            
+	            	if (dapi_roi.isEmpty())
+	            	{
+	            		System.out.println(prefix + "no nucleus");
+	            		continue;
+	            	}
+	                
+	            	
+	            	// 1.b) Improve segmentation with active contours
+	            
+	            	ActiveContours ac = new ActiveContours();
+	            	ac.input.setValue(currentField);
+	            	ac.roiInput.add(dapi_roi.toArray(new ROI[0]));
+	            	ac.edge_c.setValue(dapi_channel.getValue());
+	            	ac.edge_weight.setValue(0.2);
+	            	ac.region_c.setValue(dapi_channel.getValue());
+	            	ac.region_weight.setValue(1.0);
+	            	ac.regul_weight.setValue(0.02);
+	            	ac.division_sensitivity.setValue(0.01);
+	            	ac.evolution_bounds.setNoSequenceSelection();
+	            	ac.contour_resolution.setValue(3.0);
+	            	ac.contour_timeStep.setValue(1.0);
+	            	ac.convergence_criterion.setValue(0.3);
+	            	ac.convergence_nbIter.setValue(1000);
+	            	ac.output_roiType.setValue(ROIType.POLYGON);
+	            	ac.execute();
+	           
+	            	
+	            	// 1.c) Count the number of nuclei - proxy for the number of cells 
+	 
+	            	for (ROI roi : ac.roiOutput.getValue())
+	            	{            	            	            	
+	            		totalNuclei++;
+	            		nucleusNumber++;
+	
+	            	
+	            		// DISPLAY
+	                
+	            		((ROI2D) roi).setZ(0);
+	            		((ROI2D) roi).setT(field);
+	            		((ROI2D) roi).setC(tub_channel.getValue());
+	            		((ROI2D) roi).setColor(Color.blue);
+	            		seq.addROI(roi);
+	            	}
+		
+	            	
+	            // 2. Find MKLP1 dots (=> midbodies)
+	                  
+	            	
+	            	// 2.a) HK-Means
+	            
+	            	List<ROI> mklp1_roi = HKMeans.hKMeans(currentField, 0, mklp1_channel.getValue(), 0.0, (byte) 5,
+	            										  mklp1_minSize.getValue(), mklp1_maxSize.getValue(),
+	                    			  					  mklp1_minIntensity.getValue(), getStatus());
+	            
+	            	if (mklp1_roi.isEmpty())
+	            	{
+	            		System.out.println(prefix + "no midbody");
+	            		continue;
+	            	}
+	            
+	            	
+	            	// 2.b) Count the number of midbodies
+	                       	
+	            	totalMidbodies += mklp1_roi.size();
+	            	midbodyNumber = mklp1_roi.size();
+		              	
+	            	
+	            // 3) Count the number of cytokinesis events and remnant midbodies
+	                  
+	            	
+	            	// 3.a) Threshold tubulin signal (=> cytoplasms)
+	                    	
+	            	double[] v = KMeans.computeKMeansThresholds(currentField, tub_channel.getValue(), tub_classes.getValue(), 256);
+	            	
+	            	ROI2D tub_stain = (ROI2D) Thresholder.threshold(currentField, tub_channel.getValue(), new double[] { v[0] })[0];
+	            
+	            	
+	            		// DISPLAY            
+	            
+	            		tub_stain.setZ(0);
+	            		tub_stain.setT(field);
+	            		tub_stain.setC(tub_channel.getValue());
+	            		tub_stain.setColor(Color.green);
+	            		seq.addROI(tub_stain);
+	            
+	           
+	            	// 3.b) Classify each mklp1 ring based on its number of intersections with the surrounding stainings
+	
+	            		
+	            	// Merge the nucleic ROIs             	
+	            	
+	            	ROI2D mergeNuclei = (ROI2D) ROIUtil.merge(CollectionUtil.asList(ac.roiOutput.getValue()), BooleanOperator.OR);
+	
+	            	
+	            	// Classify each mklp1 ring 
+	            
+	            	for (int i = 0; i < mklp1_roi.size(); i++)
+	            	{
+	            		ROI2D mklp1 = (ROI2D) mklp1_roi.get(i);
+	            		
+	            		Point5D center = ROIMassCenterDescriptorsPlugin.computeMassCenter(mklp1);
+	            		double xC = MathUtil.round(center.getX(), 1);
+	            		double yC = MathUtil.round(center.getY(), 1);
+	                
+	            		Details.setRow(row, seq.getName(), field, xC, yC);
+	                
+	                
+	            		// 3.b.1) Dilate and subtract MKLP1 to get a ring
+	            		
+	            		ROI2D mklp1_outer_discus = (ROI2D) DilateROI.dilateROI(mklp1, ring_inflation.getValue(), ring_inflation.getValue(), 0);
+	            		
+	            		ROI2D mklp1_inner_discus = (ROI2D) DilateROI.dilateROI(mklp1, ring_inflation.getValue() - 2 , ring_inflation.getValue() - 2, 0);
+	                
+	            		ROI2D mklp1_ring = (ROI2D) ROIUtil.subtract(mklp1_outer_discus, mklp1_inner_discus);
+	            		
+	                
+	            		// DISPLAY 
+	                
+	            		mklp1_ring.setZ(0);
+	            		mklp1_ring.setT(field);
+	            		mklp1_ring.setC(tub_channel.getValue());
+	            		mklp1_ring.setColor(Color.red);                
+	            		seq.addROI(mklp1_ring);
+	                
+	                
+	            		// 3.b.2) Count the intersections
+	            		
+	            		int nbCrossings = 0;
+	                            		
+	            		if (!tub_stain.isEmpty() && !mklp1_ring.isEmpty() && mklp1_ring.intersects(tub_stain))
+	            		{
+	            			ROI2D intersection_tub = (ROI2D) ROIUtil.getIntersection(Arrays.asList(mklp1_ring, tub_stain));
+	            		
+	            			if (!intersection_tub.isEmpty())
+	            			{
+	            				BooleanMask2D[] possibleBridge = intersection_tub.getBooleanMask(true).getComponents();
+	                        
+	            				nbCrossings = possibleBridge.length;
+	                        
+	            			}              
+	            		}
+	                
+	            		
+	            		ROI2D intersection_dapi = (ROI2D) ROIUtil.getIntersection(Arrays.asList(mklp1_ring, mergeNuclei));
+	            		
+	            		if (!intersection_dapi.isEmpty())
+	            		{
+	            			nbCrossings = 1;
+	            		}
+	            		
+	            			
+	            		Details.setValue(row, 4, nbCrossings);
+	                
+	                
+	            		// 3.b.3) Classify each mklp1 ring
+	                
+	            		switch (nbCrossings)
+	            		{
+	                
+	            		case 0: 
+	                	
+	            			// REMNANT
+	            				
+	            			Details.setValue(row, 5, "REMNANT");
+	                    
+	            		break;
+	
+	
+	            		case 1:
+	                	
+	            			// REMNANT
+	            				
+	            			Details.setValue(row, 5, "REMNANT");
+	
+	            		break;
+	
+	
+	            		default:
+	                	
+	            			// CYTOKINESIS
+	            				
+	            			Details.setValue(row, 5, "BRIDGE");
+	                    
+	            			totalCytokinesis++;
+	            			cytokinesisNumber++;                
+	            		}
+	                 
+	            		row++;               
+	            	}
+	            
+	            FOVs.setRow(line, seq.getName(), field, nucleusNumber, midbodyNumber, cytokinesisNumber);
+	            
+	            line++;            
+	        }
+        }
+        catch (InterruptedException e)
+        {
+            // process interrupted
+        }
+	                  
+        Summary.setRow(1, seq.getName(), totalNuclei, totalMidbodies, totalCytokinesis);
+        
+        Workbooks.show(wb, "Cytokinesis");        
+    }   
+}
\ No newline at end of file
-- 
GitLab