COSEM_Predictions_Evaluation.bsh

import bdv.util.*;
import ij.*;
import ij.process.FloatProcessor;
import ij.gui.*;
import ij.Prefs;
import ij.util.ThreadUtil;
import ij.plugin.*;
import ij.process.*;
import ij.gui.*;

import java.util.concurrent.*;
import java.util.Arrays;
import java.util.ArrayList;
import java.util.stream.*;
import ij3d.*;

import org.apache.commons.math3.distribution.NormalDistribution;
import org.apache.commons.math3.distribution.BinomialDistribution;

import org.janelia.saalfeldlab.n5.*;
import org.janeliai.saalfeldlab.n5.imglib2.*;
import org.janelia.saalfeldlab.n5.imglib2.*;

import net.imglib2.img.display.imagej.*;
import net.imglib2.img.imageplus.*;
import net.imglib2.realtransform.*;
import net.imglib2.util.*;
import net.imglib2.view.*;
import net.imglib2.Cursor;


import java.awt.event.*;

import fiji.util.gui.GenericDialogPlus;

import java.io.*;
import java.nio.file.*;

//Results table
IJ.run("Set Measurements...", "area mean min redirect=None decimal=9"); //Decimal place for prediction
results_table =  new ij.measure.ResultsTable(1);
results_table.addValue("Prediction 1 Directory", "");
results_table.addValue("Prediction 2 Directory", "");
results_table.addValue("Raw Data Directory", "");
results_table.addValue("Selected Organelles", "");

//Evaluation summary plot
final Plot plot = new Plot("Evaluation Summary","Best Category","Counts");
plot.setColor("black","black");
plot.setXYLabels("{Prediction 1,Prediction 2,Neither}",null);
plot.add("separated bar",new double[]{0,0,0});
current_limits = plot.getLimits();
plot.setLimits(new double[]{current_limits[0],current_limits[1],0,1});

//Dialog boxes
data_select_dialog_box = new GenericDialogPlus("Select Data");
subvolume_setup_dialog_box = new NonBlockingGenericDialog("Subvolume Setup");
categorize_dialog_box = new NonBlockingGenericDialog("Prediction Comparison");
out_of_crops_dialog = new GenericDialogPlus("No More Crops");


//Variables
prediction_1_directory = "";
prediction_2_directory = "";
output_file = "";
training_crops_file = "";
training_crops_info = new ArrayList();

n5_prediction_1 = new N5FSReader("");
n5_prediction_2 = new N5FSReader("");
n5_raw_and_mask = new N5FSReader("");
img_mask = null;

raw_or_aligned = "";
mask_is_present = true;

start_coordinates_array = new ArrayList();
end_coordinates_array = new ArrayList();
thread_start_coordinates = new long[]{-1, -1, -1};
thread_end_coordinates = new long[]{-1, -1, -1};
start_coordinates = new long[]{-1,-1,-1};
end_coordinates = new long[]{-1,-1,-1};
crop_coordinates = "";
use_training_crops_for_evaluation = false;

organelle_checkbox_labels_list = new ArrayList();
organelle_checkbox_to_channel_mapping_list = new ArrayList();
evaluating_organelle = "";
evaluating_organelle_with_scale_1 = "";
evaluating_organelle_with_scale_2 = "";

active_channels_list = new ArrayList();
organelle_color_checkboxes = "";

raw_dimensions = new long[]{0,0,0};
raw_resolution = 0.0;
prediction_offset = new long []{0,0,0};
prediction_dimensions = new long[]{0,0,0};
prediction_1_dimensions = new long[]{0,0,0};
prediction_1_offset = new long[]{0,0,0};
prediction_2_dimensions = new long[]{0,0,0};
prediction_2_offset = new long[]{0,0,0};
prediction_resolution = 0.0;
mask_dimensions = new long[]{0,0,0};
mask_resolution = 0.0;

//scalings are with respect to highest prediction resolution
prediction_1_scaling = 1.0;
prediction_2_scaling = 1.0;
mask_scaling = 2.0;
raw_scaling = 1.0;


edge_length = 150;
organelle_image_array = new ImagePlus[2][7];
composite_images = new CompositeImage[2];

organelle_crop_count = 0;
evaluation_result_counts = new double[]{0,0,0};
evaluation_result_counts_for_stats = new double[]{0,0};

blind_image_array_assignments = new int[]{-1,-1};
composite_image_suffix = new String[]{"A","B"};

stop_coordinates_thread = false;

probability_predictions_are_this_different_by_chance = 1;


//Class for this plugin
Random thread_A_fixed_seed_generator = new Random(0); //Thread 0 fixed seed random generator
Random thread_B_fixed_seed_generator = new Random(1); //Thread 1 fixed seed random generator
Random currentThreadFixedSeedGenerator = null;
use_fixed_random_seed = false;
evaluate_segmentations = false;

public void resize(imp, width, height, depth){
		ImageProcessor processor = imp.getProcessor();
		StackProcessor sp = new StackProcessor(imp.getStack(), processor);
	    s2 = sp.resize(width, height, true);
		imp.setStack( s2);
		resizeZ(imp, depth, ImageProcessor.NONE);
}

public void resizeZ(ImagePlus imp, int newDepth, int interpolationMethod) {
        ImageStack stack1 = imp.getStack();
        int width = stack1.getWidth();
        int height = stack1.getHeight();
        int depth = stack1.getSize();
        int bitDepth = imp.getBitDepth();
        ImagePlus imp2 = IJ.createImage(imp.getTitle(), bitDepth+"-bit", width, height, newDepth);
        ImageStack stack2 = imp2.getStack();
        ImageProcessor ip = imp.getProcessor();
        ImageProcessor xzPlane1 = ip.createProcessor(width, depth);
        xzPlane1.setInterpolationMethod(0);
        ImageProcessor xzPlane2;        
        Object xzpixels1 = xzPlane1.getPixels();
        for (int y=0; y<height; y++) {
            for (int z=0; z<depth; z++) { // get xz plane at y
                Object pixels1 = stack1.getPixels(z+1);
                System.arraycopy(pixels1, y*width, xzpixels1, z*width, width);
            }
            xzPlane2 = xzPlane1.resize(width, newDepth, true);
            Object xypixels2 = xzPlane2.getPixels();
            for (int z=0; z<newDepth; z++) {
                Object pixels2 = stack2.getPixels(z+1);
                System.arraycopy(xypixels2, z*width, pixels2, y*width, width);
            }
        }
        imp.setStack(imp2.getStack());     
}
	    
public class ReadAndProcessImageInThread implements Runnable {
	//Class used to read and process images within a thread
 	
 	//Class parameters needed to read the appropriate data/prediction
 	int index_in_image_array, prediction_index;
 	double current_prediction_scaling;
 	String current_organelle;
 	n5_current = new N5FSReader("");
 	
	public ReadAndProcessImageInThread(int prediction_supplied, int index_supplied){
		//Constructor

		//Prediction and indices are provided
		prediction_index = prediction_supplied;
		index_in_image_array=index_supplied;

		//The raw data (index 0) is put in index 3 (so it is gray) in the actual image, so need to shift the lower indices
		if(index_supplied<=3){
			index_in_image_array--; 
		}
		
		if(index_supplied == 0){//raw data
			index_in_image_array=3;
			current_organelle = raw_or_aligned;
			n5_current=n5_raw_and_mask;
		}
		else{//Predictions
			current_organelle = ""; 
			if(prediction_index==0){ n5_current=n5_prediction_1; current_organelle = evaluating_organelle_with_scale_1; current_prediction_scaling = prediction_1_scaling;}
			else{n5_current=n5_prediction_2; current_organelle = evaluating_organelle_with_scale_2; current_prediction_scaling = prediction_2_scaling;}
		}
	}
	
    public void run(){
		//If an organelle is selected for this channel, read in the data, crop it and convert it to ImagePlus
		img_current = N5Utils.open(n5_current, current_organelle);
		img_current_cropped = null;
		ImagePlus imp_current_cropped  = null; 
		if (index_in_image_array==3){//This is the raw data which we adjust differently
			//might need to rescale raw data if it is at a lower resolution
			raw_start_coordinates = new long[]{Math.floor((start_coordinates[0]+prediction_offset[0])/raw_scaling), Math.floor((start_coordinates[1]+prediction_offset[1])/raw_scaling), Math.floor((start_coordinates[2]+prediction_offset[2])/raw_scaling)};
			raw_end_coordinates = new long[]{Math.floor((end_coordinates[0]+prediction_offset[0])/raw_scaling), Math.floor((end_coordinates[1]+prediction_offset[1])/raw_scaling), Math.floor((end_coordinates[2]+prediction_offset[2])/raw_scaling)};	
				
			img_current_cropped = Views.interval(img_current, raw_start_coordinates, raw_end_coordinates);				
			imp_current_cropped = ImageJFunctions.wrap(img_current_cropped, "temp");
			
			if(raw_scaling!=1){
				int width = end_coordinates[0] - start_coordinates[0] + 1;
				int height = end_coordinates[1] - start_coordinates[1] + 1;
				int depth = end_coordinates[2] - start_coordinates[2] + 1;
			
				resize(imp_current_cropped, width, height, depth);
			}
			contrastEnahncer = new ij.plugin.ContrastEnhancer();
			contrastEnahncer.stretchHistogram(imp_current_cropped, 0.35);
		}
		else{//Otherwise, we set the default to 127: need to set display range to 126 for it to start display at 127
			prediction_start_coordinates = new long[]{Math.floor(start_coordinates[0]/current_prediction_scaling), Math.floor(start_coordinates[1]/current_prediction_scaling), Math.floor(start_coordinates[2]/current_prediction_scaling)};
			prediction_end_coordinates = new long[]{Math.floor(end_coordinates[0]/current_prediction_scaling), Math.floor(end_coordinates[1]/current_prediction_scaling), Math.floor(end_coordinates[2]/current_prediction_scaling)};	
			img_current_cropped = Views.interval(img_current, prediction_start_coordinates, prediction_end_coordinates);
			imp_current_cropped = ImageJFunctions.wrap(img_current_cropped, "temp");

			if(current_prediction_scaling!=1){
				int width = end_coordinates[0] - start_coordinates[0] + 1;
				int height = end_coordinates[1] - start_coordinates[1] + 1;
				int depth = end_coordinates[2] - start_coordinates[2] + 1;
			
				resize(imp_current_cropped, width, height, depth);
			}
			
			if(evaluate_segmentations){
				imp_current_cropped.setDisplayRange(0,0.5);
			}
			else{
				imp_current_cropped.setDisplayRange(126,126);
			}
		}
		//Set the image dimensions and add it to the image array
		imp_current_cropped.setDimensions(1, edge_length, 1);
		organelle_image_array[prediction_index][index_in_image_array] = imp_current_cropped;
    }

}
 
public class MultithreadedImageReadingAndProcessing extends Thread {
	//Class to create threads and run image processing on many threads
	Thread[] threads = ThreadUtil.createThreadArray(2+2);//all organelles + 1 raw for each prediction, times 2 for two datasets
	public void run (){
		 ithread = 0;
		 for (prediction = 0; prediction<2; prediction++){//loop over the predicitions
			 for (organelle_index = 0; organelle_index<=1; organelle_index++){//loop over all the data needed for each composite; one for organelle of interest, one for predictions
			 		//Set up the actual process for each thread
					Runnable runnable = new ReadAndProcessImageInThread(prediction, organelle_index);
				   	threads[ithread] = new Thread(runnable);
				   	ithread++;
			 }
		}
		//Run the threads
		ThreadUtil.startAndJoin(threads);
	}
}


public class COSEM_Predictions_Evaluation implements PlugIn, KeyListener, MouseListener {
	public float getResolution(String path){
		String attributes =  new String(Files.readAllBytes(Paths.get(path+"/attributes.json", new String[0]))); 
		attributes = attributes.replaceAll("\\s+","");
		
        split_at = "resolution\":\\[";
        if(attributes.contains("pixelResolution")){
	   	    split_at = "\"pixelResolution\":\\{\"dimensions\":\\[";
	   	}
	   	else if(attributes.contains("transform")){
	   		split_at =  "scale\":\\[";
	   	}
	   	
		resolution_string = attributes.split(split_at)[1].split("\\]")[0].split(",");
	   	resolution = new float[]{Float.parseFloat(resolution_string[0]), Float.parseFloat(resolution_string[1]), Float.parseFloat(resolution_string[2])};
	   	return resolution[0]; //Assume symmetric resolution
	}

	public long[] getOffset(String path, float resolution){
		String attributes =  new String(Files.readAllBytes(Paths.get(path+"/attributes.json", new String[0]))); 
		attributes = attributes.replaceAll("\\s+","");
		offset = new long[]{0,0,0};
		if(attributes.contains("offset")){
        	split_at = "offset\":\\[";
        	offset_string = attributes.split(split_at)[1].split("\\]")[0].split(",");
        	offset = new long[]{Long.parseLong(offset_string[0])/resolution, Long.parseLong(offset_string[1])/resolution, Long.parseLong(offset_string[2])/resolution};
		}
		return offset;
	}
	
	//Create dialog boxes
	public void dataSelectDialogBox(){
		//Dialog box to get source directory and output file for saving results
		data_select_dialog_box.addDirectoryField("Prediction 1", Prefs.get("COSEM_Predictions_Evaluation.prediction_1_directory","/nrs/saalfeld/heinrichl/cell/gt061719/unet/02-070219/hela_cell3_314000.n5/"),75);
		data_select_dialog_box.addDirectoryField("Prediction 2", Prefs.get("COSEM_Predictions_Evaluation.prediction_2_directory","/nrs/saalfeld/heinrichl/cell/gt061719/unet/02-070219/hela_cell3_593000.n5/"),75);
		data_select_dialog_box.addDirectoryField("Raw Data",Prefs.get("COSEM_Predictions_Evaluation.raw_directory","/groups/cosem/cosem/data/HeLa_Cell3_4x4x4nm/HeLa_Cell3_4x4x4nm.n5/"),75);
		data_select_dialog_box.addFileField("Output File", Prefs.get("COSEM_Predictions_Evaluation.output_file","/groups/cosem/cosem/data/prediction_results.csv"),75);
		data_select_dialog_box.addCheckbox("Evaluate Training Crops",false);
		data_select_dialog_box.addCheckbox("Evaluate Segmentations",false);
		data_select_dialog_box.addCheckbox("Use Fixed Random Seed",false);

		data_select_dialog_box.showDialog();
		if (data_select_dialog_box.wasCanceled()) return;

		//Get directories and create corresponding n5 readers
		prediction_1_directory = data_select_dialog_box.getNextString();
		prediction_2_directory = data_select_dialog_box.getNextString();
		raw_directory = data_select_dialog_box.getNextString();
		n5_prediction_1 = new N5FSReader(prediction_1_directory);
		n5_prediction_2 = new N5FSReader(prediction_2_directory);
		n5_raw_and_mask = new N5FSReader(raw_directory);
		output_file = data_select_dialog_box.getNextString();
		use_training_crops_for_evaluation = data_select_dialog_box.getCheckboxes().get(0).getState();
		evaluate_segmentations = data_select_dialog_box.getCheckboxes().get(1).getState();
		use_fixed_random_seed = data_select_dialog_box.getCheckboxes().get(2).getState();
		if(use_fixed_random_seed) print("Using fixed random seed");
		if(evaluate_segmentations){
			print("Evaluating segmentations");
		}
		
		//Get info about training crops
		if(use_training_crops_for_evaluation){
			File[] version_directories = new File(raw_directory + "/volumes/groundtruth/").listFiles(directoryFilter);
			for(int version_index = 0; version_index<version_directories.length; version_index++){
				File[] training_crop_directories = version_directories[version_index].listFiles(directoryFilter);
				for(int crop_index=0; crop_index<training_crop_directories.length; crop_index++){
				
	    			String data =  new String(Files.readAllBytes(Paths.get(training_crop_directories[crop_index].getAbsolutePath()+"/labels/all/attributes.json", new String[0]))); 
	    			
					offset_string = data.split("\"offset\": \\[")[1].split("]")[0].replaceAll("\\s","").split(",");
					offset = new int[]{Float.parseFloat(offset_string[0])/4, Float.parseFloat(offset_string[1])/4, Float.parseFloat(offset_string[2])/4};
	
					dimensions_string = data.split("\"dimensions\": \\[")[1].split("]")[0].replaceAll("\\s","").split(",");
					dimensions = new int[]{Integer.parseInt(dimensions_string[0])/2, Integer.parseInt(dimensions_string[1])/2, Integer.parseInt(dimensions_string[2])/2};
	    			
	    			training_crops_info.add(new int[]{offset[0], offset[1], offset[2], dimensions[0], dimensions[1], dimensions[2]}); 
	    			
				}
			}
			Collections.shuffle(training_crops_info); //Randomize the order
		}
		
		
		//Add directories to results table
		results_table.addValue("Prediction 1 Directory", prediction_1_directory);
		results_table.addValue("Prediction 2 Directory", prediction_2_directory);
		results_table.addValue("Raw Data Directory", raw_directory);

		//Save to prefs
		Prefs.set("COSEM_Predictions_Evaluation.prediction_1_directory", prediction_1_directory);
		Prefs.set("COSEM_Predictions_Evaluation.prediction_2_directory", prediction_2_directory);
		Prefs.set("COSEM_Predictions_Evaluation.raw_directory", raw_directory);
		Prefs.set("COSEM_Predictions_Evaluation.output_file", output_file);
		Prefs.savePreferences();
		
		//Get image dimensions
		try{
			raw_dimensions = n5_raw_and_mask.getDatasetAttributes("/volumes/aligned").getDimensions();
			raw_or_aligned = "/volumes/aligned";
		}catch(error){
			if (new File(raw_directory+"/volumes/raw/s0").exists()){
				raw_dimensions	= n5_raw_and_mask.getDatasetAttributes("/volumes/raw/s0").getDimensions();
				raw_or_aligned = "/volumes/raw/s0";
				print("Using s0 for "+raw_directory+"/volumes/raw.");
			}
			else if(new File(raw_directory+"/volumes/raw").exists()){
				raw_dimensions	= n5_raw_and_mask.getDatasetAttributes("/volumes/raw").getDimensions();
				raw_or_aligned = "/volumes/raw";
			}
			else{
				raw_dimensions	= n5_raw_and_mask.getDatasetAttributes("/s0").getDimensions();
				raw_or_aligned = "/s0";
				print("Using s0 for "+raw_directory);
			}
		}
		raw_resolution = getResolution(raw_directory+raw_or_aligned);
		
		try{
			mask_dimensions = n5_raw_and_mask.getDatasetAttributes("/volumes/masks/foreground").getDimensions();
			//Other mask info
			img_mask = N5Utils.open(n5_raw_and_mask, "/volumes/masks/foreground");
			mask_resolution = getResolution(raw_directory+"/volumes/masks/foreground");
		}catch(error){
			print("No mask found. Proceeding without mask.");
			mask_is_present = false;
		}

	}
	
	public void subvolumeSetupDialogBox(){
		//Dialog box to select which organelles will be displayed

		//Find the organelles common to each list and create an array with them (also add None option)
		File[] predictionDirectories1 = new File(prediction_1_directory).listFiles(directoryFilter);
		File[] predictionDirectories2 = new File(prediction_2_directory).listFiles(directoryFilter);
		common_organelles_list = new ArrayList();
		for(i = 0; i<predictionDirectories1.length; i++){
			currentOrganelle = predictionDirectories1[i].getName();
			for(j=0; j<predictionDirectories2.length; j++){
				if( predictionDirectories2[j].getName().equals(currentOrganelle)){
					common_organelles_list.add(currentOrganelle);			
					break;
				}
			}
		}
		
		String[] common_organelles_array = new String[common_organelles_list.size()];
		common_organelles_array = common_organelles_list.toArray(common_organelles_array);
		subvolume_setup_dialog_box.addChoice("Organelle To Evaluate: ",common_organelles_array, common_organelles_array[0]);

		if(!use_training_crops_for_evaluation){//Then not using training crops
			//Add box for choosing image size
			subvolume_setup_dialog_box.addNumericField("Image Size (pixels)",150,0);
		}
		
		subvolume_setup_dialog_box.showDialog();

		//Add organelle choices to results table and create string for display and for table
		organelle_table_string = "";
		if (subvolume_setup_dialog_box.wasOKed()){
			organelle_choices = subvolume_setup_dialog_box.getChoices();
			evaluating_organelle = organelle_choices.get(0).getSelectedItem();
			
			if (new File(prediction_1_directory+"/"+evaluating_organelle + "/s0").exists()){
				evaluating_organelle_with_scale_1 = evaluating_organelle + "/s0";
				print("Using s0 for "+prediction_1_directory+"/"+evaluating_organelle+".");
			}
			else{
				evaluating_organelle_with_scale_1 = evaluating_organelle;
			}
			
			if (new File(prediction_2_directory+"/"+evaluating_organelle + "/s0").exists()){
				evaluating_organelle_with_scale_2 = evaluating_organelle + "/s0";
				print("Using s0 for "+prediction_2_directory+"/"+evaluating_organelle+".");
			}
			else{
				evaluating_organelle_with_scale_2 = evaluating_organelle;
			}
			
			active_channels_list.add(true);//start with that channel on
			String str = "Red: " + evaluating_organelle + " (n = 0, p-same = %.9f)   ";
			organelle_checkbox_labels_list.add(String.format(str,new Object[]{new Float(1)}));
			
			organelle_table_string += evaluating_organelle +",";
			organelle_checkbox_to_channel_mapping_list.add(0);

			organelle_checkbox_to_channel_mapping_list.add(1); //Gray will always be active by default
			organelle_checkbox_labels_list.add("Gray: Raw Data");
			organelle_table_string = organelle_table_string.substring(0,organelle_table_string.length()-1); 
			results_table.addValue("Selected Organelles", organelle_table_string);
			active_channels_list.add(true);//start with that channel on

			prediction_1_resolution = getResolution(prediction_1_directory+"/"+evaluating_organelle_with_scale_1);
			prediction_1_dimensions = n5_prediction_1.getDatasetAttributes(evaluating_organelle_with_scale_1).getDimensions(); 
			prediction_1_offset = getOffset(prediction_1_directory+"/"+evaluating_organelle_with_scale_1, prediction_1_resolution); //get offset in voxels
			prediction_2_resolution = getResolution(prediction_2_directory+"/"+evaluating_organelle_with_scale_2);
			prediction_2_dimensions = n5_prediction_2.getDatasetAttributes(evaluating_organelle_with_scale_2).getDimensions();
			prediction_2_offset = getOffset(prediction_2_directory+"/"+evaluating_organelle_with_scale_2, prediction_2_resolution); //get offset in voxels

			prediction_resolution = Math.min(prediction_1_resolution, prediction_2_resolution);
			
			prediction_1_scaling = 1.0*Math.round(1.0*prediction_1_resolution/prediction_resolution);
			prediction_2_scaling = 1.0*Math.round(1.0*prediction_2_resolution/prediction_resolution);

			// Get offset and extent in voxels by first converting to nm and then dividing by prediction_resolution which is taken to be the highest resolution
			prediction_offset = new long[]{Math.max(prediction_1_offset[0]*prediction_1_resolution, prediction_2_offset[0]*prediction_2_resolution)/prediction_resolution, 
										      Math.max(prediction_1_offset[1]*prediction_1_resolution, prediction_2_offset[1]*prediction_2_resolution)/prediction_resolution, 
										      Math.max(prediction_1_offset[2]*prediction_1_resolution, prediction_2_offset[2]*prediction_2_resolution)/prediction_resolution};
			prediction_dimensions = new long[]{Math.min(prediction_1_dimensions[0]*prediction_1_resolution, prediction_2_dimensions[0]*prediction_2_resolution)/prediction_resolution, 
											  Math.min(prediction_1_dimensions[1]*prediction_1_resolution, prediction_2_dimensions[1]*prediction_2_resolution)/prediction_resolution, 
											  Math.min(prediction_1_dimensions[2]*prediction_1_resolution, prediction_2_dimensions[2]*prediction_2_resolution)/prediction_resolution};

			//TODO: Currently we assume predictions have the same offsets and dimensions in nm
			/*prediction_dimensions = new long[]{(prediction_extents[0]-prediction_offset[0]),
											   (prediction_extents[1]-prediction_offset[1]),
											   (prediction_extents[2]-prediction_offset[2])};*/
			
			
			raw_scaling = Math.round(1.0*raw_resolution/prediction_resolution);
			mask_scaling = Math.round(1.0*mask_resolution/prediction_resolution);
			
			if(!use_training_crops_for_evaluation){
				edge_length = (Integer)subvolume_setup_dialog_box.getNextNumber();
			}

			//new Thread(test).start();
			if(!use_training_crops_for_evaluation){

				new Thread(new Runnable(){
					public void run(){
						while(!stop_coordinates_thread && start_coordinates_array.size()<=300){
							getEvaluationCropCoordinates();
						}
					}
				}).start();
				
				getValidImageCoordinates();
			}

			generateNewSubvolumes();
			results_table.show("Prediction Results");
			plot.show(); 
		}
	}
	
	FileFilter directoryFilter = new FileFilter() {
			public boolean accept(File file) {
				return file.isDirectory();
			}
	};

	public void addEvaluationToTableAndPlot(evaluation_result){
		randomness_generator = use_fixed_random_seed ? thread_A_fixed_seed_generator : ThreadLocalRandom.current();
		
		//Need to get actual prediction from blind randomization
		actual_prediction_evaluation_result = evaluation_result;
		if(evaluation_result == "Prediction A"){
			actual_prediction_evaluation_result = "Prediction  " + String.valueOf(blind_image_array_assignments[0]+1);
			evaluation_result_counts[blind_image_array_assignments[0]]+=1;
			evaluation_result_counts_for_stats[blind_image_array_assignments[0]]+=1;
		}
		else if (evaluation_result == "Prediction B"){
			actual_prediction_evaluation_result = "Prediction  " + String.valueOf(blind_image_array_assignments[1]+1);
			evaluation_result_counts[blind_image_array_assignments[1]]+=1;
			evaluation_result_counts_for_stats[blind_image_array_assignments[1]]+=1;
		}
		else{
			evaluation_result_counts[2]+=1;
			//if neither good nor bad, assign it randomly to a or b
			evaluation_result_counts_for_stats[randomness_generator.nextInt(2)]+=1;
		}

		int number_of_samples = evaluation_result_counts_for_stats[0]+evaluation_result_counts_for_stats[1];
		
		// https://stats.stackexchange.com/questions/21581/how-to-assess-whether-a-coin-tossed-900-times-and-comes-up-heads-490-times-is-bi, https://en.wikipedia.org/wiki/Binomial_distribution#Cumulative_distribution_function
		if(evaluation_result_counts_for_stats[0]==evaluation_result_counts_for_stats[1]){//otherwise, in else statement below, would double count the value
			probability_predictions_are_this_different_by_chance = 1;
		}
		else{
			binomialDistribution = new BinomialDistribution(number_of_samples,0.5);
			probability_predictions_are_this_different_by_chance = 2*binomialDistribution.cumulativeProbability((int)Math.min(evaluation_result_counts_for_stats[0], evaluation_result_counts_for_stats[1])); //times 2 since need to account for either prediction having this value
		}

		//Update plot
		plot.add("separated bar",evaluation_result_counts);
		limits = plot.getLimits();
		max_y = 0;
		for (i=0;i<2;i++){
			if (evaluation_result_counts[i]>max_y){max_y = evaluation_result_counts[i];}
		}
		new_limits = new double[]{limits[0],limits[1],limits[2],max_y+1};
		plot.setLimits(new_limits);
		plotWindow = plot.show();

		//Update results table
		results_table.setValue(evaluating_organelle+" Crop Coordinates",organelle_crop_count, crop_coordinates);  
		results_table.setValue(evaluating_organelle+" Prediction Evaluation",organelle_crop_count,actual_prediction_evaluation_result);
		results_table.setValue(evaluating_organelle+" Probability Predictions Are Same",organelle_crop_count,probability_predictions_are_this_different_by_chance);

		//Update selected organelle's checkbox label with count
		organelle_color_checkboxes = categorize_dialog_box.getCheckboxes();
		current_label = organelle_color_checkboxes.get(0).getLabel();
		old_count = String.valueOf(organelle_crop_count);
		new_count = String.valueOf(organelle_crop_count+1);
		String new_label = current_label.replace("n = " + (number_of_samples-1), "n = " +  number_of_samples );
		new_label = new_label.substring(0, new_label.length()-15) + "%.9f)   ";
		new_label = String.format(new_label, new Object[]{new Float(probability_predictions_are_this_different_by_chance)});
		organelle_color_checkboxes.get(0).setLabel(new_label);
		organelle_crop_count++;			

		//Fill rest of table and save it
		addEmptyEntriesToResultsTable(); 
		results_table.show("Prediction Results"); 
		results_table.saveAs(output_file); 

		if(!use_training_crops_for_evaluation && number_of_samples==100){
			predictions_are_different_dialog = new GenericDialogPlus("100 Comparisons Complete!");
			predictions_are_different_dialog.addMessage(String.format("After 100 comparisons, the odds that Prediction 0 and Prediction 1 would be this different by chance is: %.9f\n", new Object[]{new Float(probability_predictions_are_this_different_by_chance)}));
			predictions_are_different_dialog.showDialog();
		}

		//Get new subvolume
		generateNewSubvolumes();
	}

	public void toggleChannel(){
		//Function that toggles which channels are displayed
		for(checkbox_index=0; checkbox_index<organelle_color_checkboxes.size(); checkbox_index++){
			channel_index = organelle_checkbox_to_channel_mapping_list.get(checkbox_index);
			active_channels_list.set(channel_index, organelle_color_checkboxes.get(checkbox_index).getState());
		}
		active_channels_string = getActiveChannelsString();
		for (i=0; i<2; i++){
			composite_images[i].setActiveChannels(active_channels_string);
		}
	}

	public String getActiveChannelsString(){
		//Get currently active channels in string form required for setActiveChannels function
		active_channels_string = "";
		for (i=0; i<active_channels_list.size(); i++) {
  			channel_state = active_channels_list.get(i) ? "1" : "0";
  			active_channels_string+=channel_state;
		}
		return active_channels_string;
	}
	
	public void categorizeDialogBox(){
		//Dialog box to select which (if any) of the predictions was best
		prediction_A_button = new Button("Prediction A");
		prediction_B_button = new Button("Prediction B");
		neither_good_button = new Button("Neither - Good");
		neither_bad_button = new Button("Neither - Bad");

		categorize_dialog_box.addMessage("Evaluating Organelle: "+evaluating_organelle);

		//Add checkboxes for each organelle to turn it on/off and display evaluation counts
		String[] organelle_checkbox_labels_array = new String[organelle_checkbox_labels_list.size()];
		organelle_checkbox_labels_array = organelle_checkbox_labels_list.toArray(organelle_checkbox_labels_array);
		for(i=0; i<organelle_checkbox_labels_array.length; i++){
			categorize_dialog_box.addCheckbox(organelle_checkbox_labels_array[i],true);
		}

		//Setup callbacks to toggle channels when checkboxes toggled
		organelle_color_checkboxes = categorize_dialog_box.getCheckboxes();
		for(i=0; i<organelle_color_checkboxes.size(); i++){
			organelle_color_checkboxes.get(i).addItemListener(new ItemListener(){
				public void itemStateChanged(ItemEvent e){toggleChannel();}});
		}
		
		//Add listeners to each button to update/save the results table and get the next images
		prediction_A_button.addActionListener(new ActionListener(){ actionPerformed(ActionEvent e) { addEvaluationToTableAndPlot("Prediction A");} }); 
		prediction_B_button.addActionListener(new ActionListener(){ actionPerformed(ActionEvent e) { addEvaluationToTableAndPlot("Prediction B");} });
		neither_good_button.addActionListener(new ActionListener(){ actionPerformed(ActionEvent e) { addEvaluationToTableAndPlot("Neither - Good");} });
		neither_bad_button.addActionListener(new ActionListener(){ actionPerformed(ActionEvent e) { addEvaluationToTableAndPlot("Neither - Bad");} });

		//Add buttons and formatting
		categorize_dialog_box.addMessage("\n");//neeed so Done button actually appears
		categorize_dialog_box.add(prediction_A_button);
		categorize_dialog_box.add(prediction_B_button);
		categorize_dialog_box.add(neither_good_button);
		categorize_dialog_box.add(neither_bad_button);
		categorize_dialog_box.addMessage("\n ");//neeed so Done button actually appears
		categorize_dialog_box.setOKLabel("Done");
		categorize_dialog_box.hideCancelButton();
		categorize_dialog_box.showDialog();
	}

	public void generateNewSubvolumes(){
		//Function to actually create the merged images from a random starting coordinate. Categorize box is disabled until new images are loaded.
		
		categorize_dialog_box.disable();
		//Get image coordinates
		getValidImageCoordinates();
		//Add new coordinates to results table
		crop_coordinates = Arrays.toString(start_coordinates) + " - " + Arrays.toString(end_coordinates);
		//Run multithreaded processing
		multi = new MultithreadedImageReadingAndProcessing();
		multi.run();
		//Randomize Prediction A and B so that the user doesn't always know which side is which
		randomness_generator = use_fixed_random_seed ? thread_A_fixed_seed_generator : ThreadLocalRandom.current();
		blind_image_array_assignments[0] = randomness_generator.nextInt(2);//image corresponding to Prediction A
		blind_image_array_assignments[1] = 1-blind_image_array_assignments[0]; //image corresponding to Prediction B
		
		//Display composite images
		for (i=0; i<2; i++){//Looping over prediction A, then prediction B
			current_image_array_index = blind_image_array_assignments[i];
			
			if (composite_images[i] == null){
				composite_images[i] = ij.plugin.RGBStackMerge.mergeChannels(organelle_image_array[current_image_array_index], false);
				composite_images[i].show();
			}
			else{
				active_channels_string = getActiveChannelsString();
				composite_images[i].setImage(ij.plugin.RGBStackMerge.mergeChannels(organelle_image_array[current_image_array_index], false));
				composite_images[i].setActiveChannels(active_channels_string);			
			}
			composite_images[i].setTitle("Prediction " + composite_image_suffix[i]);
		}
		categorize_dialog_box.enable();
	}

	public void getValidImageCoordinates(){
		//Function to ensure image coordinates are in mask and have sufficient organelle content
		if(use_training_crops_for_evaluation){
			getTrainingCropCoordinates();
		}
		else{

			while(start_coordinates_array.size() == 0 || end_coordinates_array.size() == 0){ //need to wait for thread to find coordinates
				Thread.sleep(100);
			}

			start_coordinates = start_coordinates_array.get(0);
			end_coordinates = end_coordinates_array.get(0);

			start_coordinates_array.remove(0);
			end_coordinates_array.remove(0);
		}
	}

	public void getEvaluationCropCoordinates(){

		getRandomCropCoordinates();
		if(mask_is_present){
			while( !( croppedRegionIsInMask() && organelleIsInCroppedRegion())){
				getRandomCropCoordinates();
			}
		}else{
			while( !organelleIsInCroppedRegion()){
				getRandomCropCoordinates();
			}
		}
		

		start_coordinates_array.add(thread_start_coordinates);
		end_coordinates_array.add(thread_end_coordinates);
	}

	public void getTrainingCropCoordinates(){
		if(training_crops_info.size>0){
			current_crop_info = training_crops_info.get(0);
			start_coordinates[0]=current_crop_info[0];
			start_coordinates[1]=current_crop_info[1];
			start_coordinates[2]=current_crop_info[2];
			end_coordinates[0]=start_coordinates[0]+current_crop_info[3];
			end_coordinates[1]=start_coordinates[1]+current_crop_info[4];
			end_coordinates[2]=start_coordinates[2]+current_crop_info[5];
			training_crops_info.remove(0);
		}
		else{
			out_of_crops_dialog.addMessage("Ran out of crop regions to check\n");
			out_of_crops_dialog.showDialog();
		}
	}

	public void getRandomCropCoordinates(){
		randomness_generator = use_fixed_random_seed ? thread_B_fixed_seed_generator : ThreadLocalRandom.current();
		//Function to get random coordinates
		startX = randomness_generator.nextInt((Integer)prediction_dimensions[0]-(edge_length-1));
		startY = randomness_generator.nextInt((Integer)prediction_dimensions[1]-(edge_length-1));
		startZ = randomness_generator.nextInt((Integer)prediction_dimensions[2]-(edge_length-1));

		thread_start_coordinates=new long[] {startX, startY, startZ};
		thread_end_coordinates=new long[] {thread_start_coordinates[0]+edge_length-1, thread_start_coordinates[1]+edge_length-1, thread_start_coordinates[2]+edge_length-1};
	}
	
	public boolean croppedRegionIsInMask(){
		//Function to check if at least 20% of cropped region is in cell mask

		//Crop the mask which is downscaled by a factor of 2
		
		mask_start_coordinates = new long[]{Math.floor((thread_start_coordinates[0]+prediction_offset[0])/mask_scaling), Math.floor((thread_start_coordinates[1]+prediction_offset[1])/mask_scaling), Math.floor((thread_start_coordinates[2]+prediction_offset[2])/mask_scaling)};
		mask_end_coordinates = new long[]{Math.floor((thread_end_coordinates[0]+prediction_offset[0])/mask_scaling), Math.floor((thread_end_coordinates[1]+prediction_offset[1])/mask_scaling), Math.floor((thread_end_coordinates[2]+prediction_offset[2])/mask_scaling)};	
		mask_cropped = Views.interval(img_mask, mask_start_coordinates, mask_end_coordinates);
		
		ImagePlus imp_mask_cropped = ImageJFunctions.wrap(mask_cropped, "temp");
		statistics = new ij.process.StackStatistics(imp_mask_cropped,256,0,255);
		return statistics.mean>=255.0*0.2; //20% within cell
	}

	public boolean organelleIsInCroppedRegion(){
		//Function to make sure at least 10 voxels of organelle are in at least one of the subvolumes 
		cropped_region_contains_organelle = false;
		if (numValidOrganelleVoxelsInCrop(n5_prediction_1, evaluating_organelle_with_scale_1, prediction_1_scaling)>=500 || numValidOrganelleVoxelsInCrop(n5_prediction_2, evaluating_organelle_with_scale_2, prediction_2_scaling)>=500){
			cropped_region_contains_organelle = true;
		}
		return cropped_region_contains_organelle; //any voxels assigned to desired organelle
	}

	public int numValidOrganelleVoxelsInCrop(n5, evaluating_organelle_with_scale, current_prediction_scaling){
		//Function to count how many organelle voxels are in crop
		organelle = N5Utils.open(n5, evaluating_organelle_with_scale);

		prediction_start_coordinates = new long[]{Math.floor(thread_start_coordinates[0]/current_prediction_scaling), Math.floor(thread_start_coordinates[1]/current_prediction_scaling), Math.floor(thread_start_coordinates[2]/current_prediction_scaling)};
		prediction_end_coordinates = new long[]{Math.floor(thread_end_coordinates[0]/current_prediction_scaling), Math.floor(thread_end_coordinates[1]/current_prediction_scaling), Math.floor(thread_end_coordinates[2]/current_prediction_scaling)};	
		
		img_organelle_cropped = Views.interval(organelle, prediction_start_coordinates, prediction_end_coordinates);
		imp_organelle_cropped = ImageJFunctions.wrap(img_organelle_cropped, "temp");

		if(current_prediction_scaling!=1){
			int width = thread_end_coordinates[0] - thread_start_coordinates[0] + 1;
			int height = thread_end_coordinates[1] - thread_start_coordinates[1] + 1;
			int depth = thread_end_coordinates[2] - thread_start_coordinates[2] + 1;
		
			resize(imp_organelle_cropped, width, height, depth);
		}

		statistics = new ij.process.StackStatistics(imp_organelle_cropped,256,0,255);
		histogram = statistics.getHistogram();
		int num_valid_organelle_voxels = Arrays.stream(histogram, evaluate_segmentations ? 1 : 127, 256).sum(); //from is inclusive, to is exclusive so this will give us 127-255	
			
		return num_valid_organelle_voxels;
	}

	public void addEmptyEntriesToResultsTable(){
		columns = results_table.getHeadings();
		num_rows = results_table.size();
		for(row = 0; row<num_rows; row++){
			for(column = 0; column<columns.length; column++){
				if(results_table.getValueAsDouble(column, row) == 0){
					results_table.setValue(columns[column], row, "");
				}
			}
		}
	}

	public void run(){
		//Star the script
		dataSelectDialogBox();
		subvolumeSetupDialogBox();
		IJ.run("Synchronize Windows", "");
		IJ.run("Brightness/Contrast...");
		categorizeDialogBox();
		stop_coordinates_thread = true;
		WindowManager.closeAllWindows();
	}
 
}

new COSEM_Predictions_Evaluation().run(); //Run the script