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");        
    }   
}