folded_.bsh

#@ File (label="Select input directory", description="Folder with input images", style="directory") inputDir
#@ File (label="Select output directory", description="Folder to store results", style="directory") outputDir
#@ Integer(label = "Rescaling size", description="The image size will be rescaled to this value in its larger dimension", value=1000) rescaleBarImageTo
#@ Integer(label = "Scale min. pixel value", description="Minimum pixel value of pixels belonging to the scale bar", value=70) scaleMinValue
#@ Integer(label = "Scale max. pixel value", description="Maximum pixel value of pixels belonging to the scale bar", value=72) scaleMaxValue
#@ Integer(label = "Enamel min. pixel value", description="Minimum pixel value of pixels belonging to the enamel region", value=0) enamelMinValue
#@ Integer(label = "Enamel max. pixel value", description="Maximum pixel value of pixels belonging to the enamel region", value=20) enamelMaxValue
#@ Integer(label = "Dentine min. pixel value", description="Minimum pixel value of pixels belonging to the dentine region", value=210) dentineMinValue
#@ Integer(label = "Dentine max. pixel value", description="Maximum pixel value of pixels belonging to the dentine region", value=225) dentineMaxValue
#@ Integer(label = "Remaining tooth min. pixel value", description="Minimum pixel value of pixels belonging to the remaining tooth region", value=120) remToothMinValue
#@ Integer(label = "Remaining tooth max. pixel value", description="Maximum pixel value of pixels belonging to the remaining tooth region", value=140) remToothMaxValue
#@ Integer(label = "Minimum branch length", description="Minimum skeleton branch length (in pixels) to be measured", value=6) minBranchLength
#@ Boolean (label="Store images in 8-bit", description="Store when possible output images in 8-bit (32-bit otherwise)", value = true) store8bit
#@ Boolean (label="Verbose", description="Display processing information in the log window", value = false) verbose

/*
 * Script: Folded (a toolkit to describe mammalian herbivore dentition from 2D images)
 * Authors: Ignacio Arganda-Carreras (ignacio.arganda@ehu.eus), Oscar Sanisidro (oscar.sanisidro@uah.es)
 * Description:
 * 	The current script goes through all the images in the input folder and performs the following operations:
 * 		1. Rescale the image to 1000px in its larger dimension.
 * 		2. Detect, measure and remove scale bar.
 * 		3. Binarize the enamel regions based on user input values.
 * 		4. Execute LocalThickness on binarized enamel regions.
 * 		5. Execute OrientationJ Analysis on binarized enamel regions.
 * 		6. Skeletonize binarized enamel regions and extract branches.
 * 			6.1 Remove short branches if stated by user
 * 		7. Dilate the skeleton and use it as mask for the coherency and orientation images from OrientationJ.
 * 		8. Binarize the dentine and remaining tooth regions using user input values.
 * 		9. Execute OrientationJ K-means (with K = 8)
 * 		10. Execute "Fractal Box Count..." on binarized enamel region image.
 * 		11. Create an image output folder inside the output folder introduced by the user.
 * 		11. Store all resulting images and measured values (as CSV files) in the new output folder.
 */

import ij.IJ;
import ij.ImagePlus;
import ij.process.ImageConverter;
import inra.ijpb.binary.BinaryImages;
import inra.ijpb.label.LabelImages;
import inra.ijpb.segment.Threshold;
import sc.fiji.localThickness.LocalThicknessWrapper;
import ij.plugin.ImageCalculator;
import sc.fiji.analyzeSkeleton.AnalyzeSkeleton_;
import sc.fiji.analyzeSkeleton.Edge;
import sc.fiji.analyzeSkeleton.Vertex;
import sc.fiji.analyzeSkeleton.Point;
import ij.measure.ResultsTable;
import java.util.regex.Matcher;
import java.util.regex.Pattern;
import inra.ijpb.morphology.strel.DiskStrel;
import inra.ijpb.morphology.Morphology;
import inra.ijpb.measure.region2d.MaxFeretDiameter;
import inra.ijpb.measure.region2d.IntrinsicVolumesAnalyzer2D;
import inra.ijpb.algo.DefaultAlgoListener;
import orientation.imageware.Builder;
import orientation.OrientationKMeans;
import inra.ijpb.measure.IntensityMeasures;
import inra.ijpb.measure.ResultsBuilder;
import inra.ijpb.label.select.LabelSizeFiltering;
import inra.ijpb.label.select.RelationalOperator;


// Set the correct foregound and background values (for the ROIs to be correctly selected)
IJ.setForegroundColor(255, 255, 255);
IJ.setBackgroundColor(0, 0, 0);

// valid image file extensions
exts = ".tif.jpg.png.gif.tiff.jpeg.bmp.pgm.ima.dm3.dm4";
source_dir = inputDir.getAbsolutePath();
if (null == source_dir) 
{
	IJ.error("Error: No source directory was provided.");
	exit();
}

// Check if source directory exists
if( inputDir.exists() == false )
{
	IJ.error("Error: source directory " + source_dir + " does not exist.");
	exit;
}

source_dir = source_dir.replace('\\', '/');
if (!source_dir.endsWith("/")) source_dir += "/";

pathnames = inputDir.list();

for( f = 0; f < pathnames.length; f++ )
{
	IJ.showStatus( "Analyzing image " + f + " of " + pathnames.length + "..." );
	
	proportion = 1.0 * f / pathnames.length;

	pathname = pathnames[ f ];
	idot = pathname.lastIndexOf('.');
	if (-1 == idot)
		continue;
	// check extensions
	if( exts.contains(pathname.substring(idot).toLowerCase()) == false )
		continue;
	
	if( verbose )
		IJ.log( "\nProcessing image " + pathname + " ...\n" );
	image = IJ.openImage( source_dir + pathname );
	
	outputDirName = image.getShortTitle() + "_analysis";
	
	// Create output folder to store the results of the analysis of that image
	imageOutDir = new File( outputDir.getAbsolutePath() + File.separator + outputDirName );
	if( imageOutDir.exists() == false )
		imageOutDir.mkdirs();
	
	// get width and height
	imageWidth = image.getWidth();
	imageHeight = image.getHeight();
	
	// max dimension
	maxSize = imageWidth > imageHeight ? imageWidth : imageHeight;
	
	scaledInput = image.duplicate();
	
	// rescale to a max dimension (x or y) of 1000 px (defined in rescaleBarImageTo)
	if(imageHeight == imageWidth && imageHeight > rescaleBarImageTo){ 
		scaledInput = image.resize( rescaleBarImageTo, rescaleBarImageTo, "none" );
	}	
	if(imageHeight == imageWidth && imageHeight < rescaleBarImageTo){ 
		scaledInput = image.resize( rescaleBarImageTo, rescaleBarImageTo, "none" );
	}		
	if(imageHeight < imageWidth){
		newHeight = (imageHeight * rescaleBarImageTo) / imageWidth;
		newHeight = Math.round(newHeight);
		scaledInput = image.resize( rescaleBarImageTo, newHeight, "none" );
	}
	if(imageHeight > imageWidth){ 
		newWidth = (imageWidth * rescaleBarImageTo) / imageHeight;
		newWidth = Math.round(newWidth);
		scaledInput = image.resize( newWidth, rescaleBarImageTo, "none" );
	}
	
	// Get scale value from title
	pattern = Pattern.compile("scale([0-9]+)mm");
	matcher = pattern.matcher( image.getTitle() );
	
	scaleValue = -1;
	if (matcher.find()) {
		scaleValue = Double.parseDouble( matcher.group(1) );
		if( verbose )
			IJ.log( "Scale value found in file name: " + scaleValue + " mm." );
	}
	
	scaledInput.setTitle( "Scaled-"+image.getTitle() );
	//scaledInput.show();
	
	IJ.showProgress( 0.1 * proportion );
	
	// Select scale bar by color (min and max values introduced by user in HSB space)
	// The steps are as follows:
	// 1) Extract brightness channel (only if RGB)
	// 2) Apply threshold to that channel based on min and max values
	// 3) Apply connected components to extract all objects in white
	// 4) Extract largest component (that should be the scale)
	brightnessChannel = scaledInput.getProcessor();
	if( scaledInput.isRGB() )
	{
		hsbImage = scaledInput.duplicate();
		hsbImage.setTitle( "HSB-"+ scaledInput.getShortTitle() );
		// Convert image to HSB stack
		ic = new ImageConverter( hsbImage );
		ic.convertToHSB();
		hsbImage.setSlice(3); // select brightness channel
		
		brightnessChannel = hsbImage.getImageStack().getProcessor(3);
	}

	scaleBarImage = Threshold.threshold( brightnessChannel, scaleMinValue, scaleMaxValue );
	
	allLabels = BinaryImages.componentsLabeling( scaleBarImage, 4, 16 );
	labeledScale = LabelImages.keepLargestLabel( allLabels );
	
	// Create a ROI around the scale bar
	labeledScale.invert(); // invert labeled scale so "Create Selection" works
	scale = new ImagePlus( "Scale", labeledScale );
	//IJ.setThreshold( 1, 255, "raw" );
	labeledScale.setThreshold( 0, 1 );
	IJ.run( scale, "Create Selection", "");
	scaleROI = scale.getRoi();
	//scale.show();
	
	// Calculate Feret diameter using MorphoLibJ's tools
	feretCalculator = new MaxFeretDiameter();
	labels = new int[]{255};
	scale.killRoi();
	labeledScale.invert(); // invert labeled scale again so "analyzeRegions" works
	points = feretCalculator.analyzeRegions( labeledScale, labels, image.getCalibration() );
	feret = points[0].diameter();
	
	// Save Feret value in a CSV
	feretTable = new ResultsTable();
	feretTable.incrementCounter();
	feretTable.addValue("Feret", feret );	
	feretTable.saveAs( imageOutDir.getAbsolutePath() + File.separator + image.getShortTitle()+"_scale.csv" );
	
	
	// Based on that diameter, calculate pixel size
	pixelSize = scaleValue != -1 ? (scaleValue / feret) : 1;
	
	if( verbose )
		IJ.log( "Pixel size is set to " + pixelSize + " mm/pixel." );
	
	
	// Clear scale bar (set it to the same value as the first pixel in the image) in scaled image
	cleanImage = new ImagePlus( "Clean image", brightnessChannel );
	cleanImage.setRoi( scaleROI );
	IJ.run( cleanImage, "Enlarge...", "enlarge=1");
	pixelValue = image.getPixel( 0, 0 );
	IJ.setBackgroundColor( pixelValue[0], pixelValue[1], pixelValue[2] );
	IJ.run( cleanImage, "Clear", "slice" );
	IJ.setBackgroundColor( 0, 0, 0 );
	
	IJ.showProgress( 0.2 * proportion );
	
	// Binarize and save to file
	binaryIP = Threshold.threshold( brightnessChannel, enamelMinValue, enamelMaxValue );
	binary = new ImagePlus( "Enamel binary", binaryIP );
	IJ.run( binary, "Top Hat...", "radius=2 light don't");
	IJ.saveAs( binary, "Tiff", imageOutDir.getAbsolutePath() + File.separator +image.getShortTitle()+"_binary.tif" );
	
	//(new ImagePlus( "Brightness", brightnessChannel )).show();
	//binary.show();
	
	// Raw skeletonize
	referenceSkeleton = binary.duplicate();
	IJ.run( referenceSkeleton, "Skeletonize (2D/3D)", "");
	// Remove short skeleton branches
	if( minBranchLength > 0 )
	{
		if( verbose )
			IJ.log( "Removing short branches from skeleton..." );
			
		skel = new AnalyzeSkeleton_();
		skel.setup("", referenceSkeleton);
		skelResult = skel.run(AnalyzeSkeleton_.NONE, false, false, null, true, false);
		
		outStack = referenceSkeleton.getStack();
		
		// get graphs (one per skeleton in the image)
		graph = skelResult.getGraph();
		
		// list of end-points
		endPoints = skelResult.getListOfEndPoints();
		
		numRemovedBranches = 0;
		for( i = 0 ; i < graph.length; i++ )
		{
		    listEdges = graph[i].getEdges();
		
		    // go through all branches and remove branches under threshold
		    for( Edge e : listEdges )
		    {
	            if( e.getLength() < minBranchLength )
	            {
	            	p = e.getV1().getPoints();
			        v1End = endPoints.contains( p.get(0) );
			        p2 = e.getV2().getPoints();
			        v2End = endPoints.contains( p2.get(0) );
		        
		            // if any of the vertices is end-point, remove it
	                if( v1End )
	                    outStack.setVoxel( p.get(0).x, p.get(0).y, p.get(0).z, 0 );
	                if( v2End )
	                    outStack.setVoxel( p2.get(0).x, p2.get(0).y, p2.get(0).z, 0 );
	                // then, remove all slab voxels
	                for( Point p : e.getSlabs() )
	                    outStack.setVoxel( p.x, p.y, p.z, 0 );
	                numRemovedBranches++;
		        }
		    }
		}
		if( verbose )
			IJ.log( "Number of short branches removed: " + numRemovedBranches );  
	}
	
	IJ.saveAs( referenceSkeleton, "Tiff", imageOutDir.getAbsolutePath() + File.separator +image.getShortTitle()+"_skeleton.tif" );	
	
	// Thickness
	ltw = new LocalThicknessWrapper();
	thickness = ltw.processImage( binary );
	thicknessRGB = thickness.duplicate();
	//IJ.run( thicknessRGB, "Grays", "" );
	IJ.run( thicknessRGB, "RGB Color", "" ); // save as RGB color
	IJ.saveAs( thicknessRGB, "Tiff", imageOutDir.getAbsolutePath() + File.separator +image.getShortTitle()+"_thickness.tif" );	
	
	// Coherency
	dilatedSkeleton = referenceSkeleton.duplicate();
	// Dilate twice using a disk of radius 2 (diameter 5)
	dilatedIP = Morphology.dilation( dilatedSkeleton.getProcessor(), DiskStrel.fromDiameter( 5 ) );
	dilatedIP = Morphology.dilation( dilatedIP, DiskStrel.fromDiameter( 5 ) );
	dilatedSkeleton.setProcessor( dilatedIP );
	IJ.run( dilatedSkeleton, "OrientationJ Analysis", "log=0.0 tensor=6.0 gradient=4 harris-index=on orientation=on coherency=on color-survey=on s-distribution=on hue=Orientation sat=Original-Image bri=Original-Image ");
	
	// Check if OrientationJ results start with "OJ-" or not.
	prefix = IJ.getImage().getTitle().startsWith( "OJ-" ) ? "OJ-": "";
	
	IJ.selectWindow( prefix + "Orientation-1" ); 
	orientation = IJ.getImage();
	orientation.getWindow().setVisible( false );
	IJ.saveAs( orientation, "Jpeg", imageOutDir.getAbsolutePath() + File.separator +image.getShortTitle()+"_orientation.jpg" );
	
	IJ.selectWindow( prefix + "Coherency-1"); 
	coherency = IJ.getImage();
	coherency.getWindow().setVisible( false );
	
	IJ.run( coherency, "Invert", "" );
	
	dilatedSkel01 = dilatedSkeleton.duplicate(); // Make skeleton only with 0-1 values
	IJ.run( dilatedSkel01, "Divide...", "value=255");
	coherencySkel = ImageCalculator.run( dilatedSkel01, coherency, "Multiply create 32-bit" );
	if( store8bit )
	{
		grayImage = coherencySkel.duplicate();
		IJ.run( grayImage, "8-bit", "" );
		IJ.saveAs( grayImage, "Tiff", imageOutDir.getAbsolutePath() + File.separator +image.getShortTitle()+"_coherency.tif" );
	}
	else
		IJ.saveAs( coherencySkel, "Tiff", imageOutDir.getAbsolutePath() + File.separator +image.getShortTitle()+"_coherency.tif" );
	
	// Orientation
	IJ.selectWindow(prefix + "Color-survey-1"); 
	orientationRGB = IJ.getImage();
	IJ.saveAs( orientationRGB, "Tiff", imageOutDir.getAbsolutePath() + File.separator +image.getShortTitle()+"_orientationRGB.tif" );
	orientationRGB.close();
	
	//orientationSkel = ImageCalculator.run( dilatedSkel01, orientation, "Multiply create 32-bit" );
	//IJ.saveAs( orientationSkel, "Tiff", imageOutDir.getAbsolutePath() + File.separator +image.getShortTitle()+"_orientationGrayscale.tif" );
	
	IJ.showProgress( 0.3 * proportion );
	
	// K-means orientation
	orientationIW = Builder.create( orientation );
	kmeans = new OrientationKMeans();
	if( verbose )
		IJ.log( "Running orientation K-means..." );
	out = kmeans.run( orientationIW, 8, 1000 );
	kmeans = new ImagePlus("OJ-KMeans 8", out.buildImageStack() );
	IJ.selectWindow( "OJ-Table-Vector-Field-" );
	IJ.run( "Close" );
	
	// Count number of patches of each cluster
	kmeans8bit = kmeans.duplicate();
	IJ.run( kmeans8bit, "8-bit", "" ); // convert K-means image to 8 bit to use it as labels
	IJ.run( kmeans8bit, "Add...", "value=1"); // start labels in 1 and not 0
	LabelImages.remapLabels( kmeans8bit ); // remap labels to 1-8
	
	// Use tooth mask on K-means image
	toothMinValue = Math.min( Math.min( dentineMinValue, enamelMinValue ), remToothMinValue );
	toothMaxValue = Math.max( Math.max( dentineMaxValue, enamelMaxValue ), remToothMaxValue );
	toothMask = Threshold.threshold( brightnessChannel, toothMinValue, toothMaxValue );
	toothMaskIP = new ImagePlus( "Tooth mask", toothMask );
	
	IJ.run( toothMaskIP, "Divide...", "value=255");
	maskedKmeans8bit = ImageCalculator.run( toothMaskIP, kmeans8bit, "Multiply create" );
	
	// Calculate mean value of each cluster
	im = new IntensityMeasures( kmeans, maskedKmeans8bit );
	rb = new ResultsBuilder();
	rb.addResult( im.getMean() );
	
	// Calculate area of each cluster
	algo = new IntrinsicVolumesAnalyzer2D();
	algo.setDirectionNumber(4);
	algo.setConnectivity(4);
	DefaultAlgoListener.monitor(algo);
	labels8 = new int[]{1,2,3,4,5,6,7,8};
	intrinsicVolumes = algo.analyzeRegions( maskedKmeans8bit.getProcessor(), labels8, image.getCalibration() );
	
	table = new ResultsTable();
	// Go through all 8 clusters, label its patch and keep its number in a table
	lsf = new LabelSizeFiltering( RelationalOperator.GE, 11 );
	for( label = 1; label <= 8; label ++ )
	{
		table.incrementCounter();
		labelToKeep = new int[]{label};
		cluster = LabelImages.keepLabels( kmeans8bit, labelToKeep );
		clusterPatches = BinaryImages.componentsLabeling( cluster.getProcessor(), 8, 16 );
		// remove small regions of each path (<=10 pixels)
		clusterPatches = lsf.process( clusterPatches );
		LabelImages.remapLabels( clusterPatches );
		stats = clusterPatches.getStats();
		table.addValue( "Patch count", (int) stats.max );
		table.addValue( "Area (pixels^2)", intrinsicVolumes[ label-1 ].area );
	}
	rb.addResult( table );
	table = rb.getResultsTable();
	table.saveAs( imageOutDir.getAbsolutePath() + File.separator +image.getShortTitle()+"_kmeans.csv");
	
	kmeansValues = table.getColumn( "Patch count" );
	sumKmeans = 0;
	for( i=0; i < kmeansValues.length; i++ )
		sumKmeans += kmeansValues[i];
	
	IJ.showProgress( 0.7 * proportion );
	
	// Save K-means image
	IJ.run( kmeans, "physics", "" );	
	IJ.run( kmeans, "RGB Color", "" ); // save as RGB color
	IJ.run( toothMaskIP, "Multiply...", "value=255");
	IJ.run( toothMaskIP, "RGB Color", "" );
	kmeans = ImageCalculator.run( toothMaskIP, kmeans, "And create" );
	IJ.saveAs( kmeans, "Tiff", imageOutDir.getAbsolutePath() + File.separator +image.getShortTitle()+"_kmeans.tif" );
	
	// Read skeleton edge pixel values from orientation image
	
	// analyze skeleton
	skel = new AnalyzeSkeleton_();
	skel.setup("", referenceSkeleton);
	skelResult = skel.run(AnalyzeSkeleton_.NONE, false, false, null, true, false);
	 
	// get calibration
	pixelWidth = referenceSkeleton.getCalibration().pixelWidth;
	pixelHeight = referenceSkeleton.getCalibration().pixelHeight;
	pixelDepth = referenceSkeleton.getCalibration().pixelDepth;
	
	// get graphs (one per skeleton in the image)
	graph = skelResult.getGraph();
	
	if( null == graph )
	{
		IJ.error( "No enamel skeleton found!" );
		exit();
	}
	// go through all skeletons
	table = new ResultsTable();
	totalBranchLength = 0;
	totalAccumulatedDistance = 0;
	for( i = 0 ; i < graph.length; i++ )
	{
		if( verbose )
			IJ.log( " *** Skeleton " + (i+1) +" ***" );
	    listEdges = graph[i].getEdges();
	
	 	if( listEdges.size() == 0 )
	 	{
	 		if( verbose )
	 			IJ.log( "No branches found" );
	 		vertices = graph[i].getVertices();
	 		for( Vertex v : vertices )
	 		{
	 			p1 = v.getPoints().get( 0 );
	    		if( verbose ) 
	    			IJ.log( "  Point : "
		    			+ ( p1.x * pixelDepth ) + ", "
		  				+ ( p1.y * pixelHeight ) + ", "
		  				+ ( p1.z * pixelDepth ) + ", "
		  				+ orientation.getImageStack().getVoxel( p1.x, p1.y, p1.z ) );	
	 		}
	 	}
	    // go through all branches and display
	    // position in the log window
	    j=0;
	    for( Edge e : listEdges )
	    {
	    	totalBranchLength += e.getLength(); // sum lenghts to calculate OEL
	
	    	if( verbose )
	    		IJ.log( "Branch " + j + ":" );
	    	p1 = e.getV1().getPoints().get( 0 );
	    	if( verbose )
	    		IJ.log( "  Initial point : "
	    			+ ( p1.x * pixelDepth ) + ", "
	  				+ ( p1.y * pixelHeight ) + ", "
	  				+ ( p1.z * pixelDepth ) + ", "
	  				+ orientation.getImageStack().getVoxel( p1.x, p1.y, p1.z ) );				
	        
	        pList = e.getSlabs();
	        accumulatedDistanceByROI = 0;
	        for( pCount = 0; pCount < pList.size(); pCount ++  )
	        {
	        	p = pList.get( pCount );
	        	if( pCount > 0 )
	        	{
	        		dist = skel.calculateDistance( p, pList.get( pCount - 1 ) ) * pixelSize;
	        		accumulatedDistanceByROI += dist;
	        		totalAccumulatedDistance += dist;
	        	}
	        	if( verbose ) 
		  			IJ.log("  " +
		  				( p.x * pixelDepth ) + ", " +
		  				( p.y * pixelHeight ) + ", " +
		  				( p.z * pixelDepth ) + ", "
		  				+ orientation.getImageStack().getVoxel( p.x, p.y, p.z ) );	
		  		
		  		// store in table
		  		table.incrementCounter();
		  		table.addValue( "Skeleton", (i+1.0) );
		  		table.addValue( "Branch", (j+1.0) );
		  		table.addValue( "X",  p.x );
		  		table.addValue( "Y",  p.y );
		  		table.addValue( "Z",  p.z );
		  		table.addValue( "Orientation",  orientation.getImageStack().getVoxel( p.x, p.y, p.z ) );
		  		table.addValue( "Thickness (pix)",  thickness.getImageStack().getVoxel( p.x, p.y, p.z ) );
		  		table.addValue( "Thickness (mm)",  thickness.getImageStack().getVoxel( p.x, p.y, p.z ) * pixelSize );
		  		table.addValue( "Folding (1-coherency)",  coherency.getImageStack().getVoxel( p.x, p.y, p.z ) );
		  		table.addValue( "Distance by ROI (mm)",  accumulatedDistanceByROI );
		  		table.addValue( "Accumulated distance (mm)",  totalAccumulatedDistance );
	        }
	        p2 = e.getV2().getPoints().get( 0 );
	    	if( verbose )
	    		IJ.log( "  Final point : " + ( p2.x * pixelDepth ) + ", " +
	  				( p2.y * pixelHeight ) + ", " +
	  				( p2.z * pixelDepth ) + ", "
	  				+ orientation.getImageStack().getVoxel( p2.x, p2.y, p2.z ) );	
	     	j++;    
	    }
	}
	
	orientation.close();
	coherency.changes = false;
	coherency.close();
	table.saveAs( imageOutDir.getAbsolutePath() + File.separator +image.getShortTitle()+"_complexity.csv" );
	
	// Calculate mean and total folding values
	foldingValues = table.getColumn( "Folding (1-coherency)" );
	totalFolding = 0;
    for (i = 0; i < foldingValues.length; i++) {
        totalFolding += foldingValues[i];
    }
    meanFolding = totalFolding / foldingValues.length;

	// Calculate mean thickness in pixels and mm
	thicknessValues = table.getColumn( "Thickness (pix)" );
	sum = 0;
    for (i = 0; i < thicknessValues.length; i++) {
        sum += thicknessValues[i];
    }
    meanThicknessPix = sum / thicknessValues.length;
    meanThicknessMM = meanThicknessPix * pixelSize;
    
	IJ.showProgress( 0.9 * proportion );
	
	// Fractal dimension
	IJ.run( binary, "Fractal Box Count...", "box=2,3,4,6,8,12,16,32,64 black" );
	IJ.saveAs("Results", imageOutDir.getAbsolutePath() + File.separator +image.getShortTitle()+"_FractalDimension.csv");
	IJ.selectWindow( "Results" );IJ.run("Close");
	IJ.selectWindow( "Plot" );IJ.run("Close");
	
	// Caculate area and perimeter of each region
	algo = new IntrinsicVolumesAnalyzer2D();
	algo.setDirectionNumber(4);
	algo.setConnectivity(4);
	DefaultAlgoListener.monitor(algo);
	            
	// enamel region
	labeledEnamel = BinaryImages.componentsLabeling( binary.getProcessor(), 4, 16 );
	labeledEnamel = LabelImages.areaOpening( labeledEnamel, 50 ); // remove spurious regions
	labeledEnamel = Threshold.threshold( labeledEnamel, 1, 65535 );
	//(new ImagePlus("Enamel", labeledEnamel ) ).show();
	
	intrinsicVolumes = algo.analyzeRegions( labeledEnamel, labels, image.getCalibration() ); 
	table = new ResultsTable();
	enamelArea = 0;
	for( k = 0; k<intrinsicVolumes.length; k++ )
	{
		enamelArea += intrinsicVolumes[k].area;
		table.incrementCounter();
		table.addValue( "Area (pixel^2)", intrinsicVolumes[k].area );
		table.addValue( "Perimeter (pixel)", intrinsicVolumes[k].perimeter );
		table.addValue( "Area (mm^2)", intrinsicVolumes[k].area * pixelSize * pixelSize );
		table.addValue( "Perimeter (mm)", intrinsicVolumes[k].perimeter * pixelSize );
	}
	table.saveAs( imageOutDir.getAbsolutePath() + File.separator +image.getShortTitle()+"_enamel.csv");
	
	
	// dentine region
	dentineMask = Threshold.threshold( brightnessChannel, dentineMinValue, dentineMaxValue );
	labeledDentine = BinaryImages.componentsLabeling( dentineMask, 4, 16 );
	labeledDentine = LabelImages.areaOpening( labeledDentine, 50 ); // remove spurious regions
	labeledDentine = Threshold.threshold( labeledDentine, 1, 65535 );
	
	intrinsicVolumes = algo.analyzeRegions( labeledDentine, labels, image.getCalibration() ); 
	
	table = new ResultsTable();
	dentineArea = 0;
	for( k = 0; k<intrinsicVolumes.length; k++ )
	{
		dentineArea += intrinsicVolumes[k].area;
		table.incrementCounter();
		table.addValue( "Area (pixel^2)", intrinsicVolumes[k].area );
		table.addValue( "Perimeter (pixel)", intrinsicVolumes[k].perimeter );
		table.addValue( "Area (mm^2)", intrinsicVolumes[k].area * pixelSize * pixelSize );
		table.addValue( "Perimeter (mm)", intrinsicVolumes[k].perimeter * pixelSize );
	}
	table.saveAs( imageOutDir.getAbsolutePath() + File.separator +image.getShortTitle()+"_dentine.csv");
	
	// remaining tooth region
	remToothMask = Threshold.threshold( brightnessChannel, remToothMinValue, remToothMaxValue );
	labeledRemTooth = BinaryImages.componentsLabeling( remToothMask, 4, 16 );
	labeledRemTooth = LabelImages.areaOpening( labeledRemTooth, 50 ); // remove spurious regions
	labeledRemTooth = Threshold.threshold( labeledRemTooth, 1, 65535 );
	
	intrinsicVolumes = algo.analyzeRegions( labeledRemTooth, labels, image.getCalibration() ); 
	
	table = new ResultsTable();
	remainingToothArea = 0;
	for( k = 0; k<intrinsicVolumes.length; k++ )
	{
		remainingToothArea += intrinsicVolumes[k].area;
		table.incrementCounter();
		table.addValue( "Area (pixel^2)", intrinsicVolumes[k].area );
		table.addValue( "Perimeter (pixel)", intrinsicVolumes[k].perimeter );
		table.addValue( "Area (mm^2)", intrinsicVolumes[k].area * pixelSize * pixelSize );
		table.addValue( "Perimeter (mm)", intrinsicVolumes[k].perimeter * pixelSize );
	}
	table.saveAs( imageOutDir.getAbsolutePath() + File.separator +image.getShortTitle()+"_remainingTooth.csv");
	
	
	// Occlusal enamel length (OEL)
	oel = totalBranchLength;
	// Occlusal tooth area (OTA) 
	ota = enamelArea + dentineArea + remainingToothArea;
	// Occlusal Enamel Index (OEI) = OEL/√OTA
	oei = oel / Math.sqrt( ota );
	// Enamel Index (EI) EI = OEL/OTA
	ei = oel / ota;
	// Indentation index (D) D=OEL^2/(4·pi·OTA)
	d = oel * oel / ( 4 * Math.PI * ota );
	
	table = new ResultsTable();
	table.incrementCounter();
	table.addValue( "Scale bar (mm)",  scaleValue );
	table.addValue( "Scale bar (pixels)",  feret );
	table.addValue( "Occlusal enamel length (OEL) in pixels",  oel );
	table.addValue( "Occlusal tooth area (OTA) in pixels^2",  ota );
	table.addValue( "Occlusal Enamel Index (OEI)",  oei );
	table.addValue( "Enamel Index (EI)",  ei );
	table.addValue( "Indentation index (D)",  d );
	table.addValue( "2D OPC",  sumKmeans );
	
	table.addValue( "Mean folding",  meanFolding );
	table.addValue( "Total folding",  totalFolding );
	table.addValue( "Mean thickness (pix)",  meanThicknessPix );
	table.addValue( "Mean thickness (mm)",  meanThicknessMM );
	
	table.saveAs( imageOutDir.getAbsolutePath() + File.separator +image.getShortTitle()+"_summary.csv");
	
	IJ.showProgress( 1.0 * proportion );
	
}

if( verbose )
	IJ.log( "\nFinished processing of entire folder!" );