//**************************************************************************************
// Copyright (C) 2022 Mary E. Cole / Pore Extractor 2D
// First Release Date: 
//**************************************************************************************

requires("1.53g");

var filemenu = newMenu("PoreExtractor 2D Menu Tool", newArray("Image Pre-Processing", "Pore Extractor", "Pore Modifier", "Pore Analyzer", "-", "About", "Cite"));

macro "PoreExtractor 2D Menu Tool - C000D6cC000D5cD7cC000D5bC000D4aC000Da3DacC000D79C000D4bC000D7dC111D8dD93Db3C111D6bC111D6dC111Dd4C111DbcC111De5C111D8cC222D62C222D72C222D9cD9dC222D52C222D95C222D38C222C333D82C333D96C333D88C333D97C333Dc3C333D84C333D26D83C333C444D37C444D42C444D75C444Dc4De8C444D7aC444D94C444D49DdbC444D15D6aD8aD9aC444D4cDe9C555D5dC555D92C555D5aC555D39C555DdaC555D32C555D7bC555De6C555D89C555C666DcbC666DadDccC666De4C666D13C666D74D78C666De7C777D87C777D27C777D04C777D23C777C888D65C888D05C888Da2C888Df6C888D3aC888D22C888D16D99C999Dd3DeaC999D69Da6C999Db4C999D73Da5C999Da4C999Df7CaaaD63D85CaaaDaaCaaaDd5CaaaDabCaaaDf5CaaaD36CaaaD53CbbbD33D48CbbbDbbCbbbDd9CbbbDbdCbbbD3bCbbbCcccDb2CcccD43CcccD98CcccDa9CcccD7eDf8CcccD8bCcccD14CcccD6eCcccD8eCcccCdddD03CdddD4dCdddD25D66CdddD76D9bCdddDa7CdddD12D28D61CdddDdcCdddD64D9eCeeeD71CeeeD51CeeeDebCeeeD59D86CeeeDd8CeeeD81DcaCeeeD47Dc2CeeeD41CeeeD77Df4CfffDe3Df9CfffD3cD68CfffD06D5eDc5CfffDcdCfffDd6CfffD17D29D91CfffD31CfffDaeCfffD24D35Dd7CfffD02D21D46D55D56Da8Db5Dd2Dfa"{

		PE2Dmacro = getArgument();
		if (PE2Dmacro!="-") {
			if (PE2Dmacro=="Image Pre-Processing") { Preprocess(); }
			if (PE2Dmacro=="Pore Extractor") { Extract(); }
			else if (PE2Dmacro=="Pore Modifier") { Modify(); }
			else if (PE2Dmacro=="Pore Analyzer") { Analyze(); }
			else if (PE2Dmacro=="About") { About(); }
			else if (PE2Dmacro=="Cite") { Cite(); }
		}
	}

//**************************************************************************************
// Image Pre-Processing Function
// Copyright (C) 2022 Mary E. Cole / Pore Extractor 2D
// First Release Date: 
//**************************************************************************************
function Preprocess() {

setBatchMode(true);
 
requires("1.53g");

//Close any open windows or images 

Dialog.createNonBlocking("Welcome to Pore Extractor 2D!");

Dialog.addMessage("Warning: Any current images or windows will be closed! \n \nNext you will open a cross-section that is: \n \nCleared (black) within the marrow cavity\nCleared (black) outside the periosteum\n \nClick OK to begin!");

Dialog.show();

//Close existing images and windows

while (nImages>0) { 
selectImage(nImages); 
close(); 
      }

openwindows = getList("window.titles"); 
     for (i=0; i<lengthOf(openwindows); i++){ 
     window = openwindows[i]; 
     	selectWindow(window); 
     run("Close"); 
     } 


//Clear ROI manager 

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}


//Create blank variables for subfunction 

var image_choice = "";
var origimg = "";
var origname = "";
var contrast_choice = "";
var equalize_choice = "";
var elc_choice = "";
var highpass_choice = "";
var background_choice = "";
var gaussian_choice = "";
var combined_channel_choice = "";
var red_channel_choice = "";
var blue_channel_choice = "";
var green_channel_choice = "";
var dir_out = "";

//User choice of processing options 
	
Dialog.createNonBlocking("Image Preprocessing");

Dialog.setInsets(0,20,0)
Dialog.addMessage("Process Image or Folder",12,"#7f0000");

image_label = newArray("Single Image","Batch Process Folder");

Dialog.setInsets(0,25,0)
Dialog.addRadioButtonGroup("", image_label, 2, 1, "Single Image");

Dialog.setInsets(0,20,0)
Dialog.addMessage("Contrast Enhancement",12,"#7f0000");

Dialog.setInsets(0,25,0)
Dialog.addCheckbox("Enhance Global Contrast",false);
Dialog.setInsets(0,30,0)
contrast_label = newArray("Normalize","Equalize");
Dialog.addRadioButtonGroup("", contrast_label, 1, 2, "Normalize");

Dialog.setInsets(0,25,0)
Dialog.addCheckbox("Enhance Local Contrast",false);

Dialog.setInsets(0,20,0)
Dialog.addMessage("Reduce Image Noise",12,"#7f0000");

Dialog.setInsets(0,25,0)
Dialog.addCheckbox("Highpass Filter",false);

Dialog.setInsets(0,25,0)
Dialog.addCheckbox("Subtract Background",true);

Dialog.setInsets(0,25,0)
Dialog.addCheckbox("Gaussian Blur",false);

Dialog.setInsets(0,20,0)
Dialog.addMessage("Save Image Channels",12,"#7f0000");

channel_label = newArray("Combined","Red","Blue","Green");
channel_label_defaults = newArray(true,true,true,false);

Dialog.setInsets(0,25,0)
Dialog.addCheckboxGroup(4,1, channel_label, channel_label_defaults) 

Dialog.show();

image_choice = Dialog.getRadioButton();

contrast_choice = Dialog.getCheckbox();
equalize_choice = Dialog.getRadioButton();;
elc_choice = Dialog.getCheckbox();;
highpass_choice = Dialog.getCheckbox();;;
background_choice = Dialog.getCheckbox();;;;
gaussian_choice = Dialog.getCheckbox();;;;;

combined_channel_choice = Dialog.getCheckbox();;;;;;
red_channel_choice = Dialog.getCheckbox();;;;;;;
blue_channel_choice = Dialog.getCheckbox();;;;;;;;
green_channel_choice = Dialog.getCheckbox();;;;;;;;;

//Single Image----------------------------------------------------------

if(image_choice == "Single Image"){

//Prompt user to open the cleaned cross-sectional image file 

origpath = File.openDialog("Load the Cleared Brightfield Image");

//Prompt user to select output directory

dir_out=getDirectory("Select Output Location");

//Open the image and run the function 

open(origpath); 

origname=File.nameWithoutExtension;

origimg=getTitle();

Preprocess_function(origimg, origname, contrast_choice, equalize_choice, elc_choice, highpass_choice, background_choice, gaussian_choice, combined_channel_choice, red_channel_choice, blue_channel_choice, green_channel_choice, dir_out);

}

//Batch Image----------------------------------------------------------


if(image_choice == "Batch Process Folder"){

dir_in = getDirectory("Select Input Folder");

list= getFileList(dir_in); 
Array.sort(list);

//Delete all files that are not tif or bmp from file list

for (i=list.length-1; i>=0; i--) {
if (!endsWith(list[i], "tif") && !endsWith(list[i], "bmp")){
	list = Array.deleteIndex(list, i);
}
}

//If no tif or bmp files were detected, display warning

if (list.length == 0){
	exit("Macro Terminated: No .tif or .bmp images detected");
}

//Prompt user to select output directory

dir_top = getDirectory("Select Output Location");

//BEGIN IMAGE LOOP

for (i=0; i<lengthOf(list); i++) {

//Open the next image in the list 

open(dir_in+list[i]);

origname=File.nameWithoutExtension;

origimg=getTitle();

//Create the outout directory 

dir_out = dir_top+"/"+origname+"/";
File.makeDirectory(dir_out); 

Preprocess_function(origimg, origname, contrast_choice, equalize_choice, elc_choice, highpass_choice, background_choice, gaussian_choice, combined_channel_choice, red_channel_choice, blue_channel_choice, green_channel_choice, dir_out);


//End image loop 
}
//End batch process
}
//End master function call 
}


//Preprocess subfunction 


function Preprocess_function(origimg, origname, contrast_choice, equalize_choice, elc_choice, highpass_choice, background_choice, gaussian_choice, combined_channel_choice, red_channel_choice, blue_channel_choice, green_channel_choice, dir_out){

//Ensure background set to black and foreground set to white

run("Colors...", "foreground=white background=black selection=red");

//Clear any current selections on origimg 

selectImage(origimg);
run("Select None");
run("Remove Overlay");

//**CROSS-SECTIONAL GEOMETRY**

showStatus("!Extracting Borders...");

//Set measurements to area

run("Set Measurements...", "area redirect=None decimal=3");

//Threshold Total Area (TA)

selectImage(origimg);
run("Duplicate...", " ");
cortex=getTitle();

selectImage(cortex);
run("8-bit");

setAutoThreshold("Otsu dark");
setThreshold(1, 255);
setOption("BlackBackground", true);
run("Convert to Mask");

run("Clear Results");

selectImage(cortex);

run("Analyze Particles...", "display exclude clear add");

//Get maximum value of TA, in case of absolute white artifacts

roiManager("Measure");
selectWindow("Results");
setLocation(screenWidth, screenHeight);
area=Table.getColumn("Area");
Array.getStatistics(area, min, max, mean, std);
TA=(max);

var Trow = "";

run("Clear Results");
roiManager("Measure");
area=Table.getColumn("Area");

ranks = Array.rankPositions(area);
Array.reverse(ranks);

Trow=ranks[0];

//Save TA roi

roiManager("Select",Trow);
roiManager("Save", dir_out+origname+"_"+"TA.roi");
roiManager("deselect");

//Clear ROI manager 

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Clear results

run("Clear Results");

//Isolate Marrow Area (MA)

//Invert cortex by thresholding 

selectImage(cortex);
run("Select None");
run("Remove Overlay");

selectImage(cortex);
run("Invert");

run("Analyze Particles...", "display exclude clear add");

//Get maximum value of MA, in case of absolute white artifacts

roiManager("Measure");
area=Table.getColumn("Area");
Array.getStatistics(area, min, max, mean, std);
MA=(max);


var Mrow = "";

run("Clear Results");
roiManager("Measure");
area=Table.getColumn("Area");

ranks = Array.rankPositions(area);
Array.reverse(ranks);

Mrow=ranks[0];

//Save MA roi

roiManager("Select",Mrow);
roiManager("Save", dir_out+origname+"_"+"MA.roi");
roiManager("deselect");

//Clear ROI manager 

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Clear results

run("Clear Results");

selectImage(cortex);
close();

//Reopen TA and MA

roiManager("Open", dir_out+origname+"_"+"TA.roi"); 
roiManager("Open", dir_out+origname+"_"+"MA.roi"); 

//Combine TA (ROI 0) and MA (ROI 1) as CA (ROI 2)

roiManager("Select", newArray(0,1));
roiManager("XOR");

roiManager("Add");

roiManager("Select", 2);
roiManager("Save", dir_out+origname+"_"+"CA.roi");

//Remove TA and MA so that only CA remains 

roiManager("Select", newArray(0,1));
roiManager("Delete");

//Start Log 

print("**Pore Extractor 2D Preprocessing Log for " + origname + "**" );
print("");

//Enhance Global Contrast 

if(contrast_choice == 1){

showStatus("!Enhancing Global Contrast...");

if(equalize_choice == "Normalize"){
	
print("Enhance Global Contrast - Normalize Histogram");

selectImage(origimg);
roiManager("Select", 0); 

run("Enhance Contrast...", "saturated=0.3");

selectImage(origimg);
roiManager("Select", 0); 
setBackgroundColor(0, 0, 0);
run("Clear Outside");

run("Select None");
run("Remove Overlay");

}


if(equalize_choice == "Equalize"){
	
print("Enhance Global Contrast - Equalize Histogram");

selectImage(origimg);
roiManager("Select", 0); 

run("Enhance Contrast...", "saturated=0.3 equalize");

selectImage(origimg);
roiManager("Select", 0); 
setBackgroundColor(0, 0, 0);
run("Clear Outside");

run("Select None");
run("Remove Overlay");

}

}

//Enhance local contrast

if(elc_choice == 1){

print("Enhance Local Contrast");

showStatus("!Enhancing Local Contrast...");

selectImage(origimg);
run("Enhance Local Contrast (CLAHE)", "blocksize=127 histogram=256 maximum=3 mask=*None*");

selectImage(origimg);
roiManager("Select", 0); 
setBackgroundColor(0, 0, 0);
run("Clear Outside");

run("Select None");
run("Remove Overlay");

}


//Highpass filter

if(highpass_choice == 1){

print("High Pass Filter");

showStatus("!High Pass Filter...");

selectImage(origimg);
roiManager("Select", 0); 

run("Bandpass Filter...", "filter_large=40 filter_small=3 suppress=None tolerance=5 autoscale saturate");

selectImage(origimg);
roiManager("Select", 0); 
setBackgroundColor(0, 0, 0);
run("Clear Outside");

run("Select None");
run("Remove Overlay");

}

//Subtract background

if(background_choice == 1){

print("Subtract Background");

showStatus("!Subtract Background...");

selectImage(origimg);
roiManager("Select", 0); 
run("Subtract Background...", "rolling=50 light");

roiManager("Select", 0); 
setBackgroundColor(0, 0, 0);
run("Clear Outside");

run("Select None");
run("Remove Overlay");

}

//Gaussian

if(gaussian_choice == 1){

showStatus("!Gaussian Blur...");

print("Gaussian Blur");

selectImage(origimg);
run("Gaussian Blur...", "sigma=2");

selectImage(origimg);
roiManager("Select", 0); 
setBackgroundColor(0, 0, 0);
run("Clear Outside");

run("Select None");
run("Remove Overlay");

}



//If user wants no channels split or if user wants a combined image, export image as-is 

if(combined_channel_choice == 0 && red_channel_choice == 0 && blue_channel_choice == 0 && green_channel_choice == 0){

showStatus("!Saving Images...");

selectImage(origimg);
run("Select None");
run("Remove Overlay");

saveAs("TIFF", dir_out+origname+"_"+"Preprocessed.tif");
origimg=getTitle();

}


if(combined_channel_choice == 1){

showStatus("!Saving Images...");

selectImage(origimg);
run("Select None");
run("Remove Overlay");

saveAs("TIFF", dir_out+origname+"_"+"Preprocessed.tif");
origimg=getTitle();

}

//Split channels and export selected ----------------------------------------------------------------------------------------------------------------------------------------------

if(red_channel_choice == 1 || blue_channel_choice == 1 || green_channel_choice == 1){

showStatus("!Saving Images...");

print("Split Color Channels");

//Split channels 

selectImage(origimg);
run("Select None");
run("Remove Overlay");

selectImage(origimg);
run("Split Channels");

//Modify image names 

origimg_red = origimg + " (red)";
origimg_blue = origimg + " (blue)";
origimg_green = origimg + " (green)";

if(red_channel_choice == 1){

selectImage(origimg_red);
saveAs("TIFF", dir_out+origname+"_"+"Preprocessed_Red.tif");

}

if(blue_channel_choice == 1){

selectImage(origimg_blue);
saveAs("TIFF", dir_out+origname+"_"+"Preprocessed_Blue.tif");

}

if(green_channel_choice == 1){

selectImage(origimg_green);
saveAs("TIFF", dir_out+origname+"_"+"Preprocessed_Green.tif");

}

//Close bracket for split channels  
}

//Save log 

selectWindow("Log");
saveAs("Text", dir_out+origname+"_"+"Log.txt");

//Clear ROI manager 

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Clear results 

run("Clear Results");

//Close existing images and windows

while (nImages>0) { 
selectImage(nImages); 
close(); 
      }

openwindows = getList("window.titles"); 
     for (i=0; i<lengthOf(openwindows); i++){ 
     window = openwindows[i]; 
     	selectWindow(window); 
     run("Close"); 
     } 


showStatus("!Image Preprocessing Complete!");

//End preprocess function 

}


//**************************************************************************************
// Pore Extractor Function
// Copyright (C) 2022 Mary E. Cole / Pore Extractor 2D
// First Release Date: 
//**************************************************************************************

function Extract() {

setBatchMode(false);
 
requires("1.53g");

//Close existing images and windows

while (nImages>0) { 
selectImage(nImages); 
close(); 
      }

openwindows = getList("window.titles"); 
     for (i=0; i<lengthOf(openwindows); i++){ 
     window = openwindows[i]; 
     	selectWindow(window); 
     run("Close"); 
     } 
     
//Prompt user to open the cleaned cross-sectional image file 

origpath = File.openDialog("Load the Brightfield or Preprocessed Image");

//Prompt user to select output directory

dir=getDirectory("Select Output Location");
list= getFileList(dir);

//Set scale

Dialog.createNonBlocking("Set Image Scale");

Dialog.addNumber("Pixel Size:",0);
Dialog.addToSameRow();
Dialog.addMessage("pixels / mm");

Dialog.show();

scale = Dialog.getNumber();

scale_um = scale/1000;

//Open image 

open(origpath); 

origname=File.nameWithoutExtension;

origimg=getTitle();

run("Select None");
run("Remove Overlay");
setBatchMode("hide");

//Ensure background set to black and foreground set to white

run("Colors...", "foreground=white background=black selection=red");

//Clear any current  results 

run("Clear Results");


//Remove scale 

run("Set Scale...", "distance=0 known=0 pixel=1 unit=pixel");

//Start Log 

print("**Pore Extractor 2D Thresholding Log for " + origname + "**" );

print("");
print("Image Scale: " + scale + " pixels / mm");

//**CROSS-SECTIONAL GEOMETRY**

showStatus("!Extracting Borders...");

//Set measurements to area

run("Set Measurements...", "area redirect=None decimal=3");

//Threshold Total Area (TA)

selectImage(origimg);
run("Duplicate...", " ");
cortex=getTitle();
setBatchMode("hide");

selectImage(cortex);
run("8-bit");

setAutoThreshold("Otsu dark");
setThreshold(1, 255);
setOption("BlackBackground", true);
run("Convert to Mask");

run("Clear Results");

selectImage(cortex);

run("Analyze Particles...", "display exclude clear add");

//Get maximum value of TA, in case of absolute white artifacts

roiManager("Measure");
selectWindow("Results");
setLocation(screenWidth, screenHeight);
area=Table.getColumn("Area");
Array.getStatistics(area, min, max, mean, std);
TA=(max);

var Trow = "";

run("Clear Results");
roiManager("Measure");
area=Table.getColumn("Area");

ranks = Array.rankPositions(area);
Array.reverse(ranks);

Trow=ranks[0];

//Save TA roi

roiManager("Select",Trow);
roiManager("Save", dir+origname+"_"+"TA.roi");
roiManager("deselect");

//Clear ROI manager 

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Clear results

run("Clear Results");

//Isolate Marrow Area (MA)

//Invert cortex by thresholding 

selectImage(cortex);
run("Select None");
run("Remove Overlay");

selectImage(cortex);
run("Invert");

run("Analyze Particles...", "display exclude clear add");

//Get maximum value of MA, in case of absolute white artifacts

roiManager("Measure");
area=Table.getColumn("Area");
Array.getStatistics(area, min, max, mean, std);
MA=(max);


var Mrow = "";

run("Clear Results");
roiManager("Measure");
area=Table.getColumn("Area");

ranks = Array.rankPositions(area);
Array.reverse(ranks);

Mrow=ranks[0];

//Save MA roi

roiManager("Select",Mrow);
roiManager("Save", dir+origname+"_"+"MA.roi");
roiManager("deselect");

//Clear ROI manager 

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Clear results

run("Clear Results");

selectImage(cortex);
close();

//Reopen TA and MA

roiManager("Open", dir+origname+"_"+"TA.roi"); 
roiManager("Open", dir+origname+"_"+"MA.roi"); 

//Combine TA (ROI 0) and MA (ROI 1) as CA (ROI 2)

roiManager("Select", newArray(0,1));
roiManager("XOR");

roiManager("Add");

roiManager("Select", 2);
roiManager("Save", dir+origname+"_"+"CA.roi");

//Remove TA and MA so that only CA remains 

roiManager("Select", newArray(0,1));
roiManager("Delete");


//***THRESHOLDING***

Dialog.createNonBlocking("Thresholding Method");

Dialog.setInsets(5,0,0)
Dialog.addMessage("Thresholding Method",14,"#7f0000");

Dialog.setInsets(0,0,0)
thresh = newArray("Try All","Manual Pore Lumens (White)", "Manual Pore Borders (Black)", "Manual Pore Lumens + Borders", "Auto Local Phansalkar");
Dialog.addRadioButtonGroup("", thresh, 5, 1, "Try All");

Dialog.setInsets(0,13,0)
phan = 15;
Dialog.addNumber("Radius:", phan);
Dialog.addToSameRow();
Dialog.addMessage("pixels");

Dialog.show();

//Get user choices for setup

threshchoice = Dialog.getRadioButton();
phanchoice = Dialog.getNumber();

//***MANUAL GLOBAL THRESHOLDING*** - Lumens 

if (threshchoice=="Manual Pore Lumens (White)" || threshchoice=="Manual Pore Lumens + Borders" || threshchoice=="Try All"){

showStatus("!Manual Pore Lumen Thresholding...");

//Duplicate image to extract white pore contents

selectImage(origimg); 
run("Select None");
run("Remove Overlay");

selectImage(origimg); 
run("Duplicate...", "title = Manual_Pore_Lumens_(White)_Temp");
white=getTitle();

//Hide orig image

selectImage(origimg); 
setBatchMode("hide");

//Convert to 8-bit

selectImage(white);
if (bitDepth!=8)
{run("8-bit");}

//Turn on CA roi

selectImage(white);
setBatchMode("show");
roiManager("Select", 0);

//Run Auto Otsu threshold of white pore contents

setAutoThreshold("Otsu dark");
call("ij.plugin.frame.ThresholdAdjuster.setMode", "Red");
call("ij.plugin.frame.ThresholdAdjuster.setMethod","Otsu")
run("Threshold...");

//User-adjusted white pore contents thresholding

title="Pore Lumen Extraction ----------->";
msg = "Adjust the TOP slider RIGHT\n \nWhite pore lumens should be red with minimal external noise\nZoom in with the magnifying glass for precision \nThen click OK to Continue";
waitForUser(title,msg);

getThreshold(whitelower,whiteupper);
setOption("BlackBackground", true);
run("Convert to Mask");

//Clear outside black pore border image

selectImage(white);
roiManager("Select", 0); 
setBackgroundColor(0, 0, 0);
run("Clear Outside");

//Remove selection from black image

selectImage(white);
run("Select None");
run("Remove Overlay");

//Duplicate and fill holes 

selectImage(white); 
run("Duplicate...", "title = Manual_Pore_Lumens_(White)");
white_filled=getTitle(); 

//Close unfilled white image 

selectImage(white);
close();

//Fill duplicated image

selectImage(white_filled);
setBatchMode("hide");
roiManager("Select", 0); 
setOption("BlackBackground", true);
run("Close-");
run("Fill Holes");
run("Select None");
run("Remove Overlay");
selectImage(white_filled);
rename("Manual_Pore_Lumens_(White)");
white_filled=getTitle(); 

//Close threshold window

selectWindow("Threshold");
run("Close");

//Hide filled white image 

selectImage(white_filled);
setBatchMode("hide");

}


//***MANUAL GLOBAL THRESHOLDING*** - Borders

if (threshchoice=="Manual Pore Borders (Black)" || threshchoice=="Manual Pore Lumens + Borders" || threshchoice=="Try All"){


showStatus("!Manual Pore Border Thresholding...");

selectImage(origimg); 
run("Select None");
run("Remove Overlay");

selectImage(origimg); 
run("Duplicate...", "title = Manual_Pore_Borders_(Black)_Temp");
black=getTitle();

//Convert to 8-bit

selectImage(black);
if (bitDepth!=8)
{run("8-bit");}

//Turn on CA roi

selectImage(black);
setBatchMode("show");
roiManager("Select", 0);

//Run Auto Otsu threshold of black pore borders

setAutoThreshold("Otsu");
call("ij.plugin.frame.ThresholdAdjuster.setMode", "Red");
call("ij.plugin.frame.ThresholdAdjuster.setMethod","Otsu")
run("Threshold...");

//User-adjusted black pore border thresholding

title="<----------- Pore Border Extraction";
msg = "Adjust the BOTTOM slider LEFT\n \nDark pore borders should be red with minimal external noise\nZoom in with the magnifying glass for precision \nThen click OK to Continue";
waitForUser(title,msg);

//Print black thresholding to log

getThreshold(blacklower,blackupper);
setOption("BlackBackground", true);
run("Convert to Mask");

//Clear outside black pore border image

selectImage(black);
roiManager("Select", 0); 
setBackgroundColor(0, 0, 0);
run("Clear Outside");

//Remove selection from black image

selectImage(black);
run("Select None");
run("Remove Overlay");

//Duplicate and fill holes 

selectImage(black); 
run("Duplicate...", "title = Manual_Pore_Borders_(Black)");
black_filled=getTitle(); 

//Close unfilled black image 

selectImage(black);
close();

//Fill duplicated image

selectImage(black_filled);
roiManager("Select", 0); 
setOption("BlackBackground", true);
run("Close-");
run("Fill Holes");
run("Select None");
run("Remove Overlay");
selectImage(black_filled);
rename("Manual_Pore_Borders_(Black)");
black_filled=getTitle(); 

//Close threshold window

selectWindow("Threshold");
run("Close");

//Hide filled black image 

selectImage(black_filled);
setBatchMode("hide");

}


//***MANUAL GLOBAL THRESHOLDING*** - Borders + Lumens

if (threshchoice=="Manual Pore Lumens + Borders" || threshchoice=="Try All"){

	
showStatus("!Combining Images...");

//Add black pore border and white pore contents images 

imageCalculator("Add create", black_filled, white_filled);

combined=getTitle();

selectImage(combined);
run("Make Binary");

selectImage(combined);
rename("Manual_Pore_Lumens_Borders");
combined=getTitle();

selectImage(combined);
roiManager("Select", 0); 
setOption("BlackBackground", true);
run("Close-");
run("Fill Holes");
run("Select None");
run("Remove Overlay");

selectImage(combined);
setBatchMode("hide");


}


//***AUTO LOCAL THRESHOLDING***

if (threshchoice=="Auto Local Phansalkar" || threshchoice=="Try All"){

showStatus("!Auto Local Thresholding...");

selectImage(origimg);
run("Select None");
run("Remove Overlay");

selectImage(origimg); 
run("Duplicate...", "title = Auto_Local_Phansalkar");
phanimg=getTitle();

selectImage(origimg); 
setBatchMode("hide");


//Convert to 8-bit

selectImage(phanimg);
if (bitDepth!=8)
{run("8-bit");}


//Turn on CA roi

selectImage(phanimg);
roiManager("Select", 0);

//Run Phansalkar Auto Local Threshold
run("Auto Local Threshold", "method=Phansalkar radius=15 parameter_1=0 parameter_2=0");

//Clear outside auto local image

selectImage(phanimg);
roiManager("Select", 0); 
setBackgroundColor(0, 0, 0);
run("Clear Outside");

//Close and fill holes 

selectImage(phanimg);
roiManager("Select", 0); 
setOption("BlackBackground", true);
run("Close-");
run("Fill Holes");
run("Select None");
run("Remove Overlay");
selectImage(phanimg);
rename("Auto_Local_Phansalkar");
phanimg=getTitle();

}


//If user tried both - also generate a combined image of Auto Local and Combined 

if (threshchoice=="Try All"){

showStatus("!Combining Images...");

imageCalculator("Add create", combined, phanimg);

autocombined=getTitle();

selectImage(autocombined);
run("Make Binary");

selectImage(autocombined);
rename("All_Thresholds_Combined");
autocombined=getTitle();

//Close and fill holes 

selectImage(autocombined);
roiManager("Select", 0); 
setOption("BlackBackground", true);
run("Close-");
run("Fill Holes");
run("Select None");
run("Remove Overlay");

}



//**THRESHOLDED IMAGE SELECTION**

//Display options for manual global or try both 

if (threshchoice=="Manual Pore Lumens + Borders" || threshchoice=="Try All"){

//Close original image 

selectImage(origimg);
close();

//Show all hidden images 

selectImage(white_filled);
setBatchMode("show");

selectImage(black_filled);
setBatchMode("show");

selectImage(combined);
setBatchMode("show");

if (threshchoice=="Try All"){

selectImage(phanimg);
setBatchMode("show");

selectImage(autocombined);
setBatchMode("show");
	
}


//Tile windows and ask user to select best option

roiManager("Deselect"); 

run("Tile");

Dialog.createNonBlocking("Select Thresholded Image")

Dialog.addMessage("Select the threshold option with the most complete pores");

Dialog.setInsets(0,60,0)
Dialog.addImageChoice("");

Dialog.show();

origpores = Dialog.getImageChoice();

//Close all images except selected
selectImage(origpores);
close("\\Others");

//Print output to log based on user selection 


if (origpores==white_filled){
print("");
print("**Thresholding**");
print("Method: Manual Global - Otsu");
print("     White Pore Lumens: " + whitelower + " - " + whiteupper);
}


if (origpores==black_filled){
print("");
print("**Thresholding**");
print("Method: Manual Global - Otsu");
print("     Black Pore Borders: " + blacklower + " - " + blackupper);
}

if (origpores==combined){
print("");
print("**Thresholding**");
print("Combined:");
print("Method: Manual Global - Otsu");
print("     White Pore Lumens: " + whitelower + " - " + whiteupper);
print("     Black Pore Borders: " + blacklower + " - " + blackupper);
}



if (threshchoice=="Try All"){

//Print output to log based on user selection 

if (origpores==phanimg){
print("");
print("**Thresholding**");
print("Method: Auto Local - Phansalkar");
print("     Radius: " + phanchoice);
}

if (origpores==autocombined){
print("");
print("**Thresholding**");
print("Combined:");
print("Method: Manual Global - Otsu");
print("     White Pore Lumens: " + whitelower + " - " + whiteupper);
print("     Black Pore Borders: " + blacklower + " - " + blackupper);
print("Method: Auto Local - Phansalkar");
print("     Radius: " + phanchoice);
}


}

//Reopen the original image 

open(origpath); 

origname=File.nameWithoutExtension;

origimg=getTitle();

run("Select None");
run("Remove Overlay");
//setBatchMode("show");


}

//If only auto local was tried, output the log information 


if (threshchoice=="Auto Local Phansalkar"){

print("");
print("**Thresholding**");
print("Method: Auto Local - Phansalkar");
print("     Radius: " + phanchoice);

origpores = phanimg;

}

//If only pore lumens were tried, output the log information 

if (threshchoice=="Manual Pore Lumens (White)"){

print("");
print("**Thresholding**");
print("Method: Manual Global - Otsu");
print("     White Pore Lumens: " + whitelower + " - " + whiteupper);

origpores = white_filled;


}

//If only black borders were tried, output the log information 

if (threshchoice=="Manual Pore Borders (Black)"){

print("");
print("**Thresholding**");
print("Method: Manual Global - Otsu");
print("     Black Pore Borders: " + blacklower + " - " + blackupper);

origpores = black_filled;

}


//Save copy of log

selectWindow("Log");
saveAs("Text", dir+origname+"_"+"Threshold_Log.txt");
print("\\Clear");



 
//Make repeating variable 

var G_repeat = "Run";
var pore_repeat = "No, Adjust Pore Filter or Color";

//Reset origpores size
selectImage(origpores);
run("Original Scale");

//Duplicate origpores as origtemp

selectImage(origpores);
run("Duplicate...", "title=Current_Pore_Modification");
currentpores=getTitle();

selectImage(origpores);
setBatchMode("hide");

//Print Morphological modifications header


print("**Pore Extractor 2D Morphological Modification Log for " + origname + "**" );

print("");
print("Image Scale: " + scale + " pixels / mm");

print("");
print("**Morphological Modifications**");

//Set initial macro defaults; these will change to reflect whatever the person used last for Analyze Particles only 

var poresizechoice = 500;
var porecircchoice = 0.30;
var roi_color_choice = "red";
var roi_type_choice = "Outline";

var despecklechoice = false; 
var despecklecyclechoice= 1;

var outlierchoice = false;
var outliersizechoice = 1;

var closechoice = false;
var closecyclechoice = 1;

var openchoice = false;
var opencyclechoice = 1;


//Beginning of do-while loop for morphometric adjustment

do {

showStatus("!Morphological Modification...");

//Clear ROIs and overlays from previous runs

//Clear ROI manager 

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Clear overlays

selectImage(currentpores);
setBatchMode("show");
run("Select None");
run("Remove Overlay");

selectImage(origimg);
setBatchMode("show");
run("Select None");
run("Remove Overlay");

//Display Dialog
	
Dialog.createNonBlocking("Morphological Modification");

Dialog.setInsets(0,0,0)
Dialog.addMessage("Instructions",14,"#7f0000");

Dialog.setInsets(0,0,0)
Dialog.addMessage("Check one or more options to preview");

Dialog.setInsets(0,0,0)
Dialog.addMessage("You can return to this screen to select additional options");

Dialog.setInsets(0,0,0)
Dialog.addMessage("Preview Options",14,"#7f0000");

morphaction = newArray("Revert to Original Threshold Image","Preview Selected Workflow Option(s)", "Preview Current Pore ROIs");
Dialog.setInsets(0,0,0)
Dialog.addRadioButtonGroup("", morphaction, 3, 1, "Preview Current Pore ROIs");

Dialog.setInsets(0,0,0)
Dialog.addNumber("Min. Pore Size", poresizechoice);
Dialog.addToSameRow();
Dialog.addMessage(fromCharCode(956) + "m" + fromCharCode(0xb2));

Dialog.setInsets(0,0,0)
Dialog.addNumber("Min. Pore Circularity", porecircchoice);

Dialog.setInsets(0,0,0)
Dialog.addMessage("Workflow Options",14,"#7f0000");

Dialog.setInsets(0,0,0)
Dialog.addCheckbox("Despeckle", despecklechoice);
Dialog.setInsets(0,0,0)
Dialog.addToSameRow();
Dialog.addNumber("Number of Cycles", despecklecyclechoice);

Dialog.setInsets(0,0,0)
Dialog.addCheckbox("Remove Bright Outliers", outlierchoice);
Dialog.addToSameRow();
Dialog.setInsets(0,0,0)
Dialog.addNumber("Radius", outliersizechoice);
Dialog.addToSameRow();
Dialog.addMessage("pixels");

Dialog.setInsets(0,0,0)
Dialog.addCheckbox("Close and Fill Pores", closechoice);
Dialog.addToSameRow();
Dialog.setInsets(0,0,0)
Dialog.addNumber("Dilate - Erode", closecyclechoice);
Dialog.addToSameRow();
Dialog.addMessage("pixels");

Dialog.setInsets(0,0,0)
Dialog.addCheckbox("Smooth Pore Borders", openchoice);
Dialog.addToSameRow();
Dialog.setInsets(0,0,0)
Dialog.addNumber("Erode - Dilate", opencyclechoice);
Dialog.addToSameRow();
Dialog.addMessage("pixels");

Dialog.show();


previewchoice = Dialog.getRadioButton();
poresizechoice= Dialog.getNumber();
porecircchoice= Dialog.getNumber();;

despecklechoice = Dialog.getCheckbox();
despecklecyclechoice= Dialog.getNumber();;;

outlierchoice = Dialog.getCheckbox();;
outliersizechoice = Dialog.getNumber();;;;

closechoice = Dialog.getCheckbox();;;
closecyclechoice = Dialog.getNumber();;;;;

openchoice = Dialog.getCheckbox();;;;
opencyclechoice = Dialog.getNumber();;;;;;



//Radio button choice: Preview Current Pore ROIs

if (previewchoice=="Preview Current Pore ROIs"){

//Run superimpose function on unmodified current pores


do {


Dialog.createNonBlocking("Pore Filter Settings");

Dialog.setInsets(0,0,5)
Dialog.addMessage("Pore Filter Settings",14,"#7f0000");

Dialog.setInsets(0,0,0)
Dialog.addNumber("Min. Pore Size ("  + fromCharCode(956) + "m)" + fromCharCode(0xb2), poresizechoice);
//Dialog.addToSameRow();
//Dialog.addMessage(fromCharCode(956) + "m" + fromCharCode(0xb2));

Dialog.setInsets(0,0,0)
Dialog.addNumber("Min. Pore Circularity", porecircchoice);

roi_color = newArray("red", "green", "blue","magenta", "cyan", "yellow", "orange", "black", "white");

Dialog.setInsets(0,0,5)
Dialog.addMessage("Change Pore Color",14,"#7f0000");

Dialog.setInsets(0,0,0)
Dialog.addChoice("Color:", roi_color, roi_color_choice);

Dialog.setInsets(0,75,0)
Dialog.addRadioButtonGroup("",  newArray("Outline","Filled"), 1, 2, roi_type_choice);


Dialog.show();

poresizechoice= Dialog.getNumber();
porecircchoice= Dialog.getNumber();;
roi_color_choice = Dialog.getChoice();
roi_type_choice = Dialog.getRadioButton();

//Clear ROI manager 

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Clear overlays

selectImage(currentpores);
setBatchMode("show");
run("Select None");
run("Remove Overlay");

selectImage(origimg);
setBatchMode("show");
run("Select None");
run("Remove Overlay");

//Set scale 
selectImage(currentpores);
run("Set Scale...", "distance=scale_um known=1 pixel=1 unit=mm global");

//Analyze particles according to user values

selectImage(currentpores);
run("Analyze Particles...", "size=poresizechoice-Infinity circularity=porecircchoice-1.00 display exclude clear add");


//Toggle outline color

if (roi_type_choice == "Outline"){

roiManager("Deselect");
RoiManager.setPosition(0);
roiManager("Set Color", roi_color_choice);
roiManager("Set Line Width", 0);

}

//Toggle fill color

if (roi_type_choice == "Filled"){

roiManager("Deselect");
RoiManager.setPosition(0);
roiManager("Set Fill Color", roi_color_choice);

}

//Superimpose these ROIs on original cleaned image

selectImage(currentpores);
run("Scale to Fit");
selectImage(currentpores);
roiManager("Show All without labels");

selectImage(origimg);
run("Scale to Fit");
selectImage(origimg);
roiManager("Show All without labels");

roiManager("Deselect"); 

run("Tile");

setTool("zoom");

// Ask if user wants to change ROIs

Dialog.createNonBlocking("Confirm Pore Selection");

Dialog.setInsets(0,0,0)
Dialog.addRadioButtonGroup("Left-click to zoom in.\nRight-click to zoom out.\n \nAre you satisfied with final pore ROI selections?",  newArray("Yes, Save ROIs and Exit Macro","No, Adjust Pore Filter or Color","No, Adjust Morphometry"), 3, 1, "No, Adjust Pore Filter or Color");

Dialog.show();

pore_repeat = Dialog.getRadioButton();

}while (pore_repeat == "No, Adjust Pore Filter or Color");

//Once pore_repeat is to exit or adjust morphometry, proceed

if (pore_repeat == "Yes, Save ROIs and Exit Macro"){
G_repeat="Exit";
}


if (pore_repeat == "No, Adjust Morphometry"){
G_repeat="Run";	
}

//End preview current pore ROIs

}

//Radio button choice: Preview Selected Workflow Option(s)

if (previewchoice=="Preview Selected Workflow Option(s)"){

//Clear ROI manager 

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Clear overlays

selectImage(currentpores);
run("Select None");
run("Remove Overlay");

//Duplicate current pores for preview
selectImage(currentpores);
run("Duplicate...", "title=Preview_Pores");
previewpores=getTitle();

//Hide current pores 

selectImage(currentpores);
setBatchMode("hide");

//Select preview pores

selectImage(previewpores);
setBatchMode("show");

//Despeckle option 

if (despecklechoice==true){
	for(i=1; i<=despecklecyclechoice; i++){run("Despeckle");}
}

//Remove bright outliers option 

if (outlierchoice==true){
run("Remove Outliers...", "radius=outliersizechoice threshold=50 which=Bright");
}

//Morphological closing choice 

if (closechoice==true){
	run("EDM Binary Operations", "iterations=" + closecyclechoice + " operation=dilate");
	run("Close-");
	run("Fill Holes");
	run("EDM Binary Operations", "iterations=" + closecyclechoice + " operation=erode");
}

//Morphological opening choice

if (openchoice==true){
	run("EDM Binary Operations", "iterations=" + opencyclechoice + " operation=erode");
	run("EDM Binary Operations", "iterations=" + opencyclechoice + " operation=dilate");
}

//Run the ROI superimpose on the preview pores 


do {


Dialog.createNonBlocking("Pore Filter Settings");

Dialog.setInsets(0,0,5)
Dialog.addMessage("Pore Filter Settings",14,"#7f0000");

Dialog.setInsets(0,0,0)
Dialog.addNumber("Min. Pore Size ("  + fromCharCode(956) + "m)" + fromCharCode(0xb2), poresizechoice);
//Dialog.addToSameRow();
//Dialog.addMessage(fromCharCode(956) + "m" + fromCharCode(0xb2));

Dialog.setInsets(0,0,0)
Dialog.addNumber("Min. Pore Circularity", porecircchoice);

roi_color = newArray("red", "green", "blue","magenta", "cyan", "yellow", "orange", "black", "white");

Dialog.setInsets(0,0,5)
Dialog.addMessage("Change Pore Color",14,"#7f0000");

Dialog.setInsets(0,0,0)
Dialog.addChoice("Color:", roi_color, roi_color_choice);

Dialog.setInsets(0,75,0)
Dialog.addRadioButtonGroup("",  newArray("Outline","Filled"), 1, 2, roi_type_choice);


Dialog.show();

poresizechoice= Dialog.getNumber();
porecircchoice= Dialog.getNumber();;
roi_color_choice = Dialog.getChoice();
roi_type_choice = Dialog.getRadioButton();

//Clear ROI manager 

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Clear overlays

selectImage(previewpores);
run("Select None");
run("Remove Overlay");

selectImage(origimg);
setBatchMode("show");
run("Select None");
run("Remove Overlay");

//Set scale 
selectImage(previewpores);
run("Set Scale...", "distance=scale_um known=1 pixel=1 unit=mm global");

//Analyze particles according to user values

selectImage(previewpores);
run("Analyze Particles...", "size=poresizechoice-Infinity circularity=porecircchoice-1.00 display exclude clear add");


//Toggle outline color

if (roi_type_choice == "Outline"){

roiManager("Deselect");
RoiManager.setPosition(0);
roiManager("Set Color", roi_color_choice);
roiManager("Set Line Width", 0);

}

//Toggle fill color

if (roi_type_choice == "Filled"){

roiManager("Deselect");
RoiManager.setPosition(0);
roiManager("Set Fill Color", roi_color_choice);

}

//Superimpose these ROIs on original cleaned image

selectImage(previewpores);
run("Scale to Fit");
selectImage(previewpores);
roiManager("Show All without labels");

selectImage(origimg);
run("Scale to Fit");
selectImage(origimg);
roiManager("Show All without labels");

roiManager("Deselect"); 

//Hide current pores 

selectImage(currentpores);
setBatchMode("hide");

run("Tile");

setTool("zoom");

// Ask if user wants to change ROIs

Dialog.createNonBlocking("Confirm Pore Selection");

Dialog.setInsets(0,0,0)
Dialog.addRadioButtonGroup("Left-click to zoom in.\nRight-click to zoom out.\n \nAre you satisfied with final pore ROI selections?",  newArray("Yes, Save ROIs and Exit Macro","No, Adjust Pore Filter or Color","No, Adjust Morphometry"), 3, 1, "No, Adjust Pore Filter or Color");

Dialog.show();

pore_repeat = Dialog.getRadioButton();

}while (pore_repeat == "No, Adjust Pore Filter or Color");


//If user wants to exit, close currentpores and replace it with previewpores

if (pore_repeat=="Yes, Save ROIs and Exit Macro"){
selectImage(currentpores);
close();
currentpores = previewpores;

//Print changes from this workflow to log 

	if (despecklechoice==true){print("     Despeckle: " + despecklecyclechoice + " cycle(s)");}

	if (outlierchoice==true){print("     Remove Bright Outliers: " + outliersizechoice + " pixel radius");}

	if (closechoice==true){print("     Close and Fill Pores: " + closecyclechoice + " Dilate - Erode pixel(s)");}
	
	if (openchoice==true){print("     Smooth Pore Borders: " + opencyclechoice + " Erode - Dilate pixel(s)");}

G_repeat="Exit";}


//If user wants to modify morphometry, ask if they want to save the change or discard it 

if (pore_repeat=="No, Adjust Morphometry"){

Dialog.createNonBlocking("Accept Morphological Change(s)?");

Dialog.setInsets(0,0,0)
Dialog.addRadioButtonGroup("Accept Morphological Change(s)?",  newArray("Accept Change","Discard Change"), 2, 1, "Accept Change");

Dialog.show();

morphchoice = Dialog.getRadioButton();

//If user accepts change, close currentpores and replace it with previewpores

if (morphchoice=="Accept Change"){
selectImage(currentpores);
close();
currentpores = previewpores;

//Show current pores 

selectImage(currentpores);
setBatchMode("show");

//Print changes from this workflow to log 

	if (despecklechoice==true){print("     Despeckle: " + despecklecyclechoice + " cycle(s)");}

	if (outlierchoice==true){print("     Remove Bright Outliers: " + outliersizechoice + " pixel radius");}

	if (closechoice==true){print("     Close and Fill Pores: " + closecyclechoice + " Dilate - Erode pixel(s)");}
	
	if (openchoice==true){print("     Smooth Pore Borders: " + opencyclechoice + " Erode - Dilate pixel(s)");}

}

//If user discards change, close previewpores

if (morphchoice=="Discard Change"){
	
selectImage(previewpores);
close();

}

G_repeat="Run";

//End No, Adjust Morphometry Loop

}

//End Preview Selected Workflow Option(s)

}


//Radio button choice: Revert to Original Threshold Image

if (previewchoice=="Revert to Original Threshold Image"){

//Close current image and re-duplicate from original thresholded 

selectImage(currentpores);
close();
selectImage(origpores);
run("Duplicate...", "title=Current_Pore_Modification");
currentpores=getTitle();

//Restart log

print("\\Clear");

print("**Pore Extractor 2D Morphological Modification Log for " + origname + "**" );

print("");
print("Image Scale: " + scale + " pixels / mm");

print("");
print("**Morphological Modifications**");

//Trigger morphometry loop restart

G_repeat="Run";

}


//End of do-while loop for morphometric adjustment

}while (G_repeat=="Run");


//Print Analyze Particle Parameters

print("");
print("**Analyze Particle Parameters**");
print("     Minimum Pore Size: " + poresizechoice + " mm" + fromCharCode(0xb2));
print("     Minimum Pore Circularity: " + porecircchoice);


//Save ROI set with user color choice

roiManager("Deselect");
RoiManager.setPosition(0);
roiManager("Set Color", roi_color_choice);
roiManager("Set Line Width", 0);

roiManager("Save", dir+origname+"_"+"PoreExtractor_RoiSet.zip");

//Save copy of log

selectWindow("Log");
saveAs("Text", dir+origname+"_"+"Morphological_Log.txt");

//Close all windows and images

while (nImages>0) { 
selectImage(nImages); 
close(); 
      }

openwindows = getList("window.titles"); 
     for (i=0; i<lengthOf(openwindows); i++){ 
     window = openwindows[i]; 
     	selectWindow(window); 
     run("Close"); 
     } 

showStatus("!Pore Extraction Complete!");

}

//**************************************************************************************
// Pore Modifier Function
// Copyright (C) 2022 Mary E. Cole / Pore Extractor 2D
// First Release Date: 
//**************************************************************************************


function Modify() {

requires("1.53g");

setBatchMode(false);

//Close all windows and images

while (nImages>0) { 
selectImage(nImages); 
close(); 
      }

openwindows = getList("window.titles"); 
     for (i=0; i<lengthOf(openwindows); i++){ 
     window = openwindows[i]; 
     	selectWindow(window); 
     run("Close"); 
     } 

//Prompt user to open the cleaned cross-sectional image file 

origpath = File.openDialog("Load the Brightfield or Preprocessed Image");

open(origpath); 

origimg=getTitle();

//Save to IJ Preferences

origname=File.nameWithoutExtension;

call("ij.Prefs.set", "orig.string", origname);

//Clear overlays

selectImage(origimg);
run("Select None");
run("Remove Overlay");

//Remove scale 

run("Set Scale...", "distance=0 known=0 pixel=1 unit=pixel");

//Prompt user to load ROI set to modify 

roipath = File.openDialog("Open the Pore ROIset To Modify");

//Clear ROI manager 

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Open ROI Path

roiManager("open", roipath);

//Prompt user to select output directory

dir=getDirectory("Select Output Location");

//Save to IJ Preferences

call("ij.Prefs.set", "appendroi.dir", dir);

//Check whether Previous labels already exist

	nR = roiManager("Count"); 
	roiIdxorig = newArray(nR); 
	p=0; 
	 
	for (i=0; i<nR; i++) { 
		roiManager("Select", i); 
		rName = Roi.getName(); 
//Find indecies of original ROIs
		if (startsWith(rName, "Prev_") ) { 
			roiIdxorig[p] = i; 
			p++; 
		} 
	} 

//If this is not a modified ROI set, then change all  to Prev_
	
roipathname = File.getNameWithoutExtension(roipath);

if (!matches(roipathname, ".*ModifiedRoiSet.*")){

for (i = 0; i<roiManager("count"); i++){
		roiManager("select", i);
		roiManager("rename", "Prev_"+i)}

}

//Display ROIs on loaded image 

roiManager("Show All without labels");

//Set original colors for selection and save

var roi_color_choice = "red";
call("ij.Prefs.set", "roi_color_choice.string", roi_color_choice );

var roi_color_save = "green";
call("ij.Prefs.set", "roi_color_save.string", roi_color_save);

//Change selection color to original choice 

run("Colors...", "foreground=white background=black selection="+ roi_color_choice);


//Do-While loop for manual modification

do{

//Instructions popup

Dialog.createNonBlocking("Manual Modification of ROI Set");

Dialog.setInsets(0,0,0)
Dialog.addMessage("Removing Incorrect ROIs",14,"#7f0000");

Dialog.setInsets(0,0,0)
Dialog.addMessage("Toggle labels ON [F1] and click label of incorrect ROI")

Dialog.setInsets(0,0,0)
Dialog.addMessage("Press backspace key to delete")

Dialog.setInsets(5,0,0)
Dialog.addMessage("Adding Missed ROIs",14,"#7f0000");

Dialog.setInsets(0,0,0)
Dialog.addMessage("Toggle labels OFF [F2] to view selection quality")

Dialog.setInsets(0,0,0)
Dialog.addMessage("Select Wand tool [1] and modulate tolerance [2] to auto-select ROIs")

Dialog.setInsets(0,0,0)
Dialog.addMessage("Or use Freehand Selection tool [3] to outline ROIs with mouse / stylus")

Dialog.setInsets(0,0,0)
Dialog.addMessage("Press t key to add selection to ROI manager")

Dialog.setInsets(5,0,0)
Dialog.addMessage("Modifying Existing ROIs with the Selection Brush",14,"#7f0000");

Dialog.setInsets(0,0,0)
Dialog.addMessage("Selection brush diameter can be decreased [F7] or increased [F8]")

Dialog.setInsets(0,0,0)
Dialog.addMessage("Use keyboard shortcut [6] and click label of ROI to be modified")

Dialog.setInsets(0,0,0)
Dialog.addMessage("Use keyboard shortcut [6] again to switch to selection brush tool")

Dialog.setInsets(0,0,0)
Dialog.addMessage("Split ROI: ")

Dialog.setInsets(0,0,0)
Dialog.addMessage("           Hold down [Alt] and draw along split line, then use [F5] to split")

Dialog.setInsets(0,0,0)
Dialog.addMessage("Expand ROI: ")

Dialog.setInsets(0,0,0)
Dialog.addMessage("            Hold down [Shift] and draw along rim of ROI to be added, then use [F6] to fill hole")


Dialog.setInsets(0,0,0)
Dialog.addMessage("Revert ROI: ")

Dialog.setInsets(0,0,0)
Dialog.addMessage("            [Ctrl] + [Shift] + [E] without clicking outside selection")


Dialog.setInsets(5,0,0)
Dialog.addMessage("Keyboard Shortcuts",14,"#7f0000");

Dialog.setInsets(0,0,0)
Dialog.addMessage("[1]    Wand Tool")

Dialog.setInsets(0,0,0)
Dialog.addMessage("[2]    Adjust Wand Tool Tolerance")

Dialog.setInsets(0,0,0)
Dialog.addMessage("[3]    Freehand Selection Tool (Draw with Mouse/Stylus)")

Dialog.setInsets(0,0,0)
Dialog.addMessage("[4]    Zoom Tool (Left-Click = Zoom In, Right-Click = Zoom Out) ")

Dialog.setInsets(0,0,0)
Dialog.addMessage("[5]    Scrolling Tool (Grab and Drag)")

Dialog.setInsets(0,0,0)
Dialog.addMessage("[6]    Selection Brush Tool")

Dialog.setInsets(0,0,0)
Dialog.addMessage("[F1]   Toggle ROI Labels On")

Dialog.setInsets(0,0,0)
Dialog.addMessage("[F2]   Toggle ROI Labels Off")

Dialog.setInsets(0,0,0)
Dialog.addMessage("[F3]   Toggle ROIs Off")

Dialog.setInsets(0,0,0)
Dialog.addMessage("[F4]   Reset Wand Tool Tolerance to Zero")

Dialog.setInsets(0,0,0)
Dialog.addMessage("[F5]   Split ROI After Selection Brush Division")

Dialog.setInsets(0,0,0)
Dialog.addMessage("[F6]   Fill ROI After Selection Brush Expansion")

Dialog.setInsets(0,0,0)
Dialog.addMessage("[F7]   Decrease Selection Brush Size by 5 pixels")

Dialog.setInsets(0,0,0)
Dialog.addMessage("[F8]   Increase Selection Brush Size by 5 pixels")

Dialog.setInsets(0,0,0)
Dialog.addMessage("[F9]   Revert ROI Immediately After Split or Fill")

Dialog.setInsets(0,0,0)
Dialog.addMessage("[F10]  Save New Version of ROI Set")


Dialog.setInsets(5,0,0)
Dialog.addMessage("Change Pore Colors",14,"#7f0000");

roi_color = newArray("red", "green", "blue","magenta", "cyan", "yellow", "orange", "black", "white");

Dialog.setInsets(0,0,0)
Dialog.addChoice("Original Pore Color:", roi_color, roi_color_choice);

Dialog.setInsets(0,0,0)
Dialog.addChoice("Modified Pore Color:", roi_color, roi_color_save);

Dialog.setInsets(5,0,0)
Dialog.addMessage("Exiting Macro",14,"#7f0000");

Dialog.setInsets(0,0,0)
modifydialog = newArray("Click OK to update colors","Click OK to save final ROI set and exit");
Dialog.addRadioButtonGroup("", modifydialog, 2, 1, "Click OK to update colors");

Dialog.setInsets(0,0,0)
Dialog.addMessage("ROI Color Code",14,"#7f0000");

//Dark Color Options 

if (roi_color_choice == "red"){roi_color_choice_display = "#bf0000";}
if (roi_color_choice == "green"){roi_color_choice_display = "#007f04";}
if (roi_color_choice == "blue"){roi_color_choice_display = "#0a13c2";}
if (roi_color_choice == "magenta"){roi_color_choice_display = "#c20aaf";}
if (roi_color_choice == "cyan"){roi_color_choice_display = "#0ac2a0";}
if (roi_color_choice == "yellow"){roi_color_choice_display = "#ab9422";} 
if (roi_color_choice == "orange"){roi_color_choice_display = "#b87d25";}
if (roi_color_choice == "black"){roi_color_choice_display = "black";}
if (roi_color_choice == "white"){roi_color_choice_display = "#a8a6a3";}


if (roi_color_save == "red"){roi_color_save_display = "#bf0000";}
if (roi_color_save == "green"){roi_color_save_display = "#007f04";}
if (roi_color_save == "blue"){roi_color_save_display = "#0a13c2";}
if (roi_color_save == "magenta"){roi_color_save_display = "#c20aaf";}
if (roi_color_save == "cyan"){roi_color_save_display = "#0ac2a0";}
if (roi_color_save == "yellow"){roi_color_save_display = "#ab9422";} 
if (roi_color_save == "orange"){roi_color_save_display = "#b87d25";}
if (roi_color_save == "black"){roi_color_save_display = "black";}
if (roi_color_save == "white"){roi_color_save_display = "#a8a6a3";}


Dialog.setInsets(0,5,0)
Dialog.addMessage("Unmodified ROIs",12,roi_color_choice_display);
Dialog.setInsets(0,5,0)
Dialog.addMessage("Modified ROIs",12,roi_color_save_display);

Dialog.show();

exitmodify = Dialog.getRadioButton();

roi_color_choice = Dialog.getChoice();
call("ij.Prefs.set", "roi_color_choice.string", roi_color_choice );

roi_color_save = Dialog.getChoice();
call("ij.Prefs.set", "roi_color_save.string", roi_color_save);

//Change selection color to color choice 

run("Colors...", "foreground=white background=black selection="+ roi_color_choice);

//Get array of updated ROIs

//Check whether the updated ROIs include the name previous 
//Print index of ROIs that do not include the name previous 
//Modified from: https://forum.image.sc/t/selecting-roi-based-on-name/3809

	nR = roiManager("Count"); 
	roiIdx = newArray(nR); 
	roiIdxorig = newArray(nR); 
	k=0; 
	p=0;
	clippedIdx = newArray(0); 
	 
	for (i=0; i<nR; i++) { 
		roiManager("Select", i); 
		rName = Roi.getName(); 
//Find indecies of new ROIs
		if (!startsWith(rName, "Prev_") ) { 
			roiIdx[k] = i; 
			k++; 
		} 
//Find indecies of original ROIs
		if (startsWith(rName, "Prev_") ) { 
			roiIdxorig[p] = i; 
			p++; 
		} 
	} 
//Make array of new ROIs
if (k>0) { 
	prevIdx = Array.trim(roiIdxorig,p);
	clippedIdx = Array.trim(roiIdx,k); 
	} 
//Make array of new ROIs
if (k==0) { 
	prevIdx = Array.trim(roiIdxorig,p); 
	} 


//Change original and previous ROIs to user choice

for (i = 0; i<prevIdx.length; i++){
		roiManager("select",prevIdx[i]);
		roiManager("Set Color", roi_color_choice);}

for (i = 0; i<clippedIdx.length; i++){
		roiManager("select",clippedIdx[i]);
		roiManager("Set Color", roi_color_save);}

} while (exitmodify=="Click OK to update colors");

//Save final ROI and exit macro 

roiRename(dir,origname,roi_color_choice,roi_color_save);

//Delete temp file roi if it exists

temproidir = dir+origname+"_"+"Temp"+".roi";

if (File.exists(temproidir))
{File.delete(temproidir);}

//Close existing images and windows

while (nImages>0) { 
selectImage(nImages); 
close(); 
      }

openwindows = getList("window.titles"); 
     for (i=0; i<lengthOf(openwindows); i++){ 
     window = openwindows[i]; 
     	selectWindow(window); 
     run("Close"); 
     } 

showStatus("!Pore Modification Complete!");

exit();


}



//**************************************************************************************
// Keyboard Shortcuts for Pore Modifier
// Copyright (C) 2022 Mary E. Cole / Pore Extractor 2D
// First Release Date: 
//**************************************************************************************

//Keyboard shortcuts-----------------------------------------
//Numpad [n] shortcuts require NumLock


//Keyboard shortcuts-----------------------------------------
//Numpad [n] shortcuts require NumLock

macro "Wand Tool [1]" {setTool("wand");}

	macro "Wand Tool [n1]" {setTool("wand");}

macro "Wand Options [2]" {run("Wand Tool...");}

	macro "Wand Options [n2]" {run("Wand Tool...");}

macro "Freehand Selection Tool [3]" {setTool("freehand");}

	macro "Freehand Selection Tool [n3]" {setTool("freehand");}

macro "Zoom Tool [4]" {setTool("zoom");}

	macro "Zoom Tool [n4]" {setTool("zoom");}

macro "Scrolling Tool [5]" {setTool("hand");}

	macro "Scrolling Tool [n5]" {setTool("hand");}

macro "Selection Brush Tool [6]" {
//Check whether ROIs are selected
currentroi = roiManager("index");
//If no ROIs are selected, switch to freehand selection tool 
if (currentroi < 0){
	//Turn on labels
	if (roiManager("count")==0) showStatus("ROI Manager is empty");
	else (roiManager("Show All with labels"));
	setTool("freehand");}
//If ROI is selected, switch to selection brush
if (currentroi >= 0){setTool("brush");}   
}

	macro "Selection Brush Tool [n6]" {
	//Check whether ROIs are selected
	currentroi = roiManager("index");
	//If no ROIs are selected, switch to freehand selection tool 
	if (currentroi < 0){setTool("freehand");}
	//If ROI is selected, switch to selection brush
	if (currentroi >= 0){setTool("brush");}   
	}


macro "ROI Labels On [f1]" {
	if (roiManager("count")==0) showStatus("ROI Manager is empty");
	else (roiManager("Show All with labels"));
}

macro "ROI Labels Off [f2]" {
	if (roiManager("count")==0) showStatus("ROI Manager is empty");
	else (roiManager("Show All without labels"));
}

macro "ROIs Off [f3]" {
	if (roiManager("count")==0) showStatus("ROI Manager is empty");
	else (roiManager("Show None"));
}

macro "Reset Wand Tolerance [f4]" {
	setTool("wand");
	run("Wand Tool...", "tolerance=0");
}

macro "Split ROI [f5]" {

dirpath= call("ij.Prefs.get",  "appendroi.dir", "nodir");
origpath= call("ij.Prefs.get", "orig.string", "");

//Save the current ROI 

roiManager("Save", dirpath+origpath+"_"+"Temp"+".roi");

//Save the current last ROI index

roicount = roiManager("count") - 1;

call("ij.Prefs.set", "roi.number", roicount);

//Get current ROI index
currentroi = roiManager("index");

roiManager("Update");
roiManager("Split");
//Delete original ROI
roiManager("Select", currentroi);
roiManager("Delete");
run("Select None");

//Show the revised image 
roiManager("Show All");


}

macro "Fill ROI [f6]" {

dir= call("ij.Prefs.get",  "appendroi.dir", "nodir");
origname= call("ij.Prefs.get", "orig.string", "");

//Save the current ROI 

roiManager("Save", dir+origname+"_"+"Temp"+".roi");

//Save the current last ROI index

roicount = roiManager("count") - 1;

call("ij.Prefs.set", "roi.number", roicount);

//Save the current ROI 

roiManager("Save", dir+origname+"_"+"Temp"+".roi");

//Set binary options to black background

setOption("BlackBackground", true);

//Get current ROI index
currentroi = roiManager("index");

roiManager("Update");

run("Create Mask");
setBatchMode("hide")

selectWindow("Mask");

run("Fill Holes");

selectWindow("Mask");
run("Analyze Particles...", "exclude include add");

selectWindow("Mask");
close();
setBatchMode("show")

//Delete original ROI
roiManager("Select", currentroi);
roiManager("Delete");
run("Select None");

//Show the revised image 
roiManager("Show All");

}


macro "Decrease Selection Brush Size [f7]" {
 
    setTool("brush");

 	a=call("ij.gui.Toolbar.getBrushSize");
 	step_in = Math.abs(a);
 	step_out = step_in - 5;
 	//Floor is 5 pixels 
 	step_out = Math.max(step_out, 5)
    b=call("ij.gui.Toolbar.setBrushSize",step_out);
    a=call("ij.gui.Toolbar.getBrushSize");
    showStatus("Selection Brush Size: " + step_out);
    
}


macro "Increase Selection Brush Size [f8]" {
 
    setTool("brush");
 	a=call("ij.gui.Toolbar.getBrushSize");
 	step_in = Math.abs(a);
 	step_out = step_in + 5;
    b=call("ij.gui.Toolbar.setBrushSize",step_out);
    a=call("ij.gui.Toolbar.getBrushSize");
    showStatus("Selection Brush Size: " + step_out);
   
}

macro "Revert ROI [f9]" {

roiManager("Deselect");

roicurrent = call("ij.Prefs.get", "roi.number", "");

roiend = roiManager("count");

roidiff = Math.abs(roiend)-Math.abs(roicurrent);

var roibound = newArray;

for (i = 0; i<roidiff; i++){
	roibound [i] = Math.abs(roicurrent)+i;}

//Delete all new rois formed after the split or expansion

roiManager("Select", roibound);

roiManager("Delete");

//Reopen saved ROI 

dir= call("ij.Prefs.get",  "appendroi.dir", "nodir");
origname= call("ij.Prefs.get", "orig.string", "");

roiManager("Open", dir+origname+"_"+"Temp"+".roi");

//Show the revised image 
roiManager("Show All");

}

macro "Save ROIs [f10]" {
	dirpath= call("ij.Prefs.get",  "appendroi.dir", "nodir");
	origpath= call("ij.Prefs.get", "orig.string", "");
	roi_color_choice_path = call("ij.Prefs.get", "roi_color_choice.string", "");
	roi_color_save_path = call("ij.Prefs.get", "roi_color_save.string", "");
	roiRename(dirpath,origpath,roi_color_choice_path,roi_color_save_path);
  }


//Function to append ROI names-------------------------- 

function roiRename(dir,origname,roi_color_choice,roi_color_save){

//Check whether the updated ROIs include the name previous 
//Print index of ROIs that do not include the name previous 
//Modified from: https://forum.image.sc/t/selecting-roi-based-on-name/3809

	nR = roiManager("Count"); 
	roiIdx = newArray(nR); 
	roiIdxorig = newArray(nR); 
	k=0; 
	p=0;
	clippedIdx = newArray(0); 
	 
	for (i=0; i<nR; i++) { 
		roiManager("Select", i); 
		rName = Roi.getName(); 
//Find indecies of new ROIs
		if (!startsWith(rName, "Prev_") ) { 
			roiIdx[k] = i; 
			k++; 
		} 
//Find indecies of original ROIs
		if (startsWith(rName, "Prev_") ) { 
			roiIdxorig[p] = i; 
			p++; 
		} 
	} 
//Make array of new ROIs
if (k>0) { 
	prevIdx = Array.trim(roiIdxorig,p);
	clippedIdx = Array.trim(roiIdx,k); 
	} 
//Make array of new ROIs
if (k==0) { 
	prevIdx = Array.trim(roiIdxorig,p); 
	} 

//Change original and previous ROIs to user choice

for (i = 0; i<prevIdx.length; i++){
		roiManager("select",prevIdx[i]);
		roiManager("Set Color", roi_color_choice);}

for (i = 0; i<clippedIdx.length; i++){
		roiManager("select",clippedIdx[i]);
		roiManager("Set Color", roi_color_save);}
		
//Pad function from https://imagej.nih.gov/ij/macros/misc/Conference%20Macros/07_Functions.ijm

function leftPad(n, width) {
      s =""+n;
      while (lengthOf(s)<width)
          s = "0"+s;
      return s;
  }

//Get current file list from directory

list1 = getFileList(dir);
Array.sort(list1);


for (j=list1.length-1; j>=0; j--) {
if (!matches(list1[j], ".*ModifiedRoiSet.*")){
	list1 = Array.deleteIndex(list1, j);
}
}


//If there are no existing ModifiedRoiSet files, start with 001

if (lengthOf(list1) == 0){
	roiappend = 1;
	roiappendpaste = leftPad(roiappend,3);
}

//If there are existing ModifiedRoiSet files, find the highest value

if (lengthOf(list1) > 0){

Array.sort(list1);

roilast = list1[lengthOf(list1)-1];

splitend = split(roilast, "_");

roiend = splitend[lengthOf(splitend)-1];

splitfile = split(roiend, ".");

roinum = splitfile[0];

roinumint = parseInt(roinum);

roiappend = roinumint+1;

roiappendpaste = leftPad(roiappend,3);

}

//Save the roi set with the new append number 

roiManager("Deselect");
roiManager("Save", dir+origname+"_"+"ModifiedRoiSet"+"_"+roiappendpaste+".zip");

//Re-display ROIs

roiManager("Show All with labels")


}




//-----------------------------------------------------------------------------------------------------------------------------------------


//**************************************************************************************
// Pore Analyzer Function
// Copyright (C) 2022 Mary E. Cole / Pore Extractor 2D
// First Release Date: 
//**************************************************************************************


function Analyze() {

requires("1.53g");

setBatchMode(true);

setBatchMode("hide");

//Close existing images and windows

while (nImages>0) { 
selectImage(nImages); 
close(); 
      }

openwindows = getList("window.titles"); 
     for (i=0; i<lengthOf(openwindows); i++){ 
     window = openwindows[i]; 
     	selectWindow(window); 
     run("Close"); 
     } 

//Clear any current results 

run("Clear Results");


//Check for BoneJ Slice Geometry Installation 

List.setCommands;
    
if (List.get("Slice Geometry")=="") {
    Dialog.createNonBlocking("BoneJ Slice Geometry Not Detected");
    Dialog.addMessage("Use Help --> Update --> Manage Update Sites --> BoneJ --> Close --> Apply Changes");
    Dialog.addMessage("Then restart ImageJ and retry");
	Dialog.show();
	exit("BoneJ Installation Required");
	}

     
//Prompt user to open the cleaned cross-sectional image file 

origpath = File.openDialog("Load the Brightfield or Preprocessed Image");

//Prompt user to load ROI set to modify 

roipath = File.openDialog("Open Final Pore ROIset");

//Prompt user to select output directory

dir=getDirectory("Select Output Location");

//Set scale

Dialog.createNonBlocking("Set Image Scale");

Dialog.addNumber("Pixel Size:",0);
Dialog.addToSameRow();
Dialog.addMessage("pixels / mm");

Dialog.show();

scale = Dialog.getNumber();
scale_um = scale/1000;

//Open the cleaned cross-sectional image 

open(origpath); 

origname=File.nameWithoutExtension;

origimg=getTitle();

//Remove any overlays

selectImage(origimg);
run("Select None");
run("Remove Overlay");

//Ensure background set to black and foreground set to white

run("Colors...", "foreground=white background=black selection=red");

//Remove scale 

run("Set Scale...", "distance=0 known=0 pixel=1 unit=pixel");

//***CROSS-SECTIONAL BOUNDING***********************************************************************************************************
//Modified from Pore Extractor Function 
//Will save over TA, MA, and CA images and ROIs if they exist in the output folder in case the user has modified the original image

showStatus("!Cross-Sectional Geometry...");

//Create folder subdivision for cross-sectional output

geomdir=dir+"/Cross-Sectional Geometry/";
File.makeDirectory(geomdir); 

//Set measurements to area

run("Set Measurements...", "area redirect=None decimal=3");

//Threshold Total Area (TA)

selectImage(origimg);
run("Duplicate...", " ");
cortex=getTitle();

selectImage(cortex);
run("8-bit");

setAutoThreshold("Otsu dark");
setThreshold(1, 255);
setOption("BlackBackground", true);
run("Convert to Mask");

run("Clear Results");

selectImage(cortex);

run("Analyze Particles...", "display exclude clear add");

//Get maximum value of TA, in case of absolute white artifacts

roiManager("Measure");
selectWindow("Results");
setLocation(screenWidth, screenHeight);
        
area=Table.getColumn("Area");
Array.getStatistics(area, min, max, mean, std);
TA=(max);

var Trow = "";

run("Clear Results");
roiManager("Measure");
area=Table.getColumn("Area");

ranks = Array.rankPositions(area);
Array.reverse(ranks);

Trow=ranks[0];

//Save TA roi

roiManager("Select",Trow);
roiManager("Save", geomdir+origname+"_"+"TA.roi");
roiManager("deselect");

//Make blank image of TA and save in output directory

selectImage(origimg);
run("Duplicate...", " ");
total=getTitle();

//Clear total slice

selectImage(total);
run("Select All");
setBackgroundColor(0, 0, 0);
run("Clear", "slice");

//Select and flatten TA image

selectImage(total);
roiManager("Select", Trow);
run("Fill", "slice");

//Clear ROI manager 

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}


selectImage(total);
run("Select None");
run("Remove Overlay");

selectImage(cortex);
run("Select None");
run("Remove Overlay");

//Threshold and save image 

selectImage(total);
run("8-bit");
setAutoThreshold("Otsu");
//run("Threshold...");
setThreshold(1, 255);
setOption("BlackBackground", true);
run("Convert to Mask");

saveAs("Tiff", geomdir+origname+"_"+"TA.tif");
total=getTitle();


//Clear results

run("Clear Results");

//Isolate Marrow Area (MA)

//Invert cortex by thresholding 

selectImage(cortex);
run("8-bit");

setAutoThreshold("Otsu");
//run("Threshold...");
setThreshold(1, 255);
setOption("BlackBackground", true);
run("Convert to Mask");

run("Invert");

run("Analyze Particles...", "display exclude clear add");

//Get maximum value of MA, in case of absolute white artifacts

roiManager("Measure");
area=Table.getColumn("Area");
Array.getStatistics(area, min, max, mean, std);
MA=(max);


var Mrow = "";

run("Clear Results");
roiManager("Measure");
area=Table.getColumn("Area");

ranks = Array.rankPositions(area);
Array.reverse(ranks);

Mrow=ranks[0];

//Save MA roi

roiManager("Select",Mrow);
roiManager("Save", geomdir+origname+"_"+"MA.roi");
roiManager("deselect");

//Clear outside marrow 

selectImage(cortex);
roiManager("Select", Mrow);
run("Fill", "slice");
run("Clear Outside");

//Clear ROI manager 

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

selectImage(cortex);
run("Select None");
run("Remove Overlay");

//Clear results

run("Clear Results");

//Threshold and save image 

selectImage(cortex);
run("8-bit");
setAutoThreshold("Otsu");
//run("Threshold...");
setThreshold(1, 255);
setOption("BlackBackground", true);
run("Convert to Mask");

selectImage(cortex);
saveAs("Tiff",geomdir+origname+"_"+"MA.tif");
marrow=getTitle();


//Combine MA and TA to get CA image

selectImage(marrow);

imageCalculator("Subtract create",total,marrow);

cortical=getTitle();

run("8-bit");
setAutoThreshold("Otsu");
//run("Threshold...");
setThreshold(1, 255);
setOption("BlackBackground", true);
run("Convert to Mask");

saveAs("Tiff", geomdir+origname+"_"+"CA.tif");
cortical=getTitle();

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Clear results

run("Clear Results");

//Reopen TA and MA

roiManager("Open", geomdir+origname+"_"+"TA.roi"); 
roiManager("Open", geomdir+origname+"_"+"MA.roi"); 

//Combine TA (ROI 0) and MA (ROI 1) as CA (ROI 2)

roiManager("Select", newArray(0,1));
roiManager("XOR");

roiManager("Add");

roiManager("Select", 2);
roiManager("Save", geomdir+origname+"_"+"CA.roi");

//BONEJ ANALYSIS ON CA

//Clear results 

run("Clear Results");

//Set scale for cross-sectional measurements according to user input

run("Set Scale...", "distance=scale known=1 pixel=1 unit=mm global");

//Run BoneJ Slice Geometry

selectImage(cortical);
run("Select None");
run("Remove Overlay");

selectImage(cortical);
run("Slice Geometry", "bone=unknown bone_min=1 bone_max=255 slope=0.0000 y_intercept=0");

//Get row numbers of columns

selectWindow("Results");
text = getInfo();
lines = split(text, "\n");
columns = split(lines[0], "\t");

if (columns[0]==" ")
 columns[0]= "Number";

//Pull variables from BoneJ results table

 Imin= getResult(columns[13],0);
 Imax = getResult(columns[14],0);

 Zpol = getResult(columns[18],0);

//**CROSS-SECTIONAL MORPHOMETRY***********************************************************************************************************

//Make table for exporting RCA and PI values

rca="[Relative Cortical Area]"; 
run("New... ", "name="+rca+" type=Table"); 
print(rca,"\\Headings:Image\tTotal Area (mm^2)\tMarrow Area (mm^2)\tCortical Area (mm^2)\t% Cortical Area\t% Bone Area Total Porosity Corrected\t% Bone Area Cortical Porosity Corrected\t% Bone Area Trabecularized Porosity Corrected\t% Marrow Area\tParabolic Index (Y)\tImin (mm^4)\tImax (mm^4)\tZpol (mm^3)\t");
selectWindow("Relative Cortical Area");
setLocation(screenWidth, screenHeight);

//Set scale for cross-sectional measurements according to user input

run("Set Scale...", "distance=scale known=1 pixel=1 unit=mm global");

//Set measurements to area only 

run("Set Measurements...", "area redirect=None decimal=3");

run("Clear Results");

//Measure TA, MA, and CA

roiManager("Select", 0);
roiManager("Measure");
TA=getResult("Area");

roiManager("Select", 1);
roiManager("Measure");
MA=getResult("Area");

roiManager("Select", 2);
roiManager("Measure");
CA=getResult("Area");

//Clear results 

run("Clear Results");

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Compute RCA and PI

RCA=(CA/TA)*100;

PerMA = (MA/TA)*100;

Para=(CA*MA)/(TA*TA);

//Clear results 

run("Clear Results");

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}



//**DERIVE TOTAL PORE IMAGE**

//Create folder subdivision for whole cross-sectional output 

poredir=dir+"/Pore Types/";
File.makeDirectory(poredir); 

//Derive pore image from ROIs and re-derive because flattening causes edge chatter

//Make blank duplicate to hold output total pore image

selectImage(marrow);  
run("Select None");
run("Remove Overlay");

selectImage(marrow); 
run("Duplicate...", "title=[TotTemp]");
pores="TotTemp";

selectWindow(pores);
run("Select All");
setBackgroundColor(0, 0, 0);
run("Clear", "slice");
run("Clear Outside", "slice");
run("Select None");
run("Remove Overlay");

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}


//Open final pore ROI set via user-provided path

roiManager("open", roipath);

//Renumber

for (i = 0; i<roiManager("count"); i++){
		roiManager("select", i);
		roiManager("rename", i)}

//Flatten total pores and save

	roiManager("deselect");
	roiManager("fill")
	selectWindow(pores);
	roiManager("Show All without Labels");
	run("Flatten");
	selectWindow("TotTemp-1");
	rename("TotalPores"); 
	poresfin="TotalPores";
	
selectWindow(poresfin);
run("8-bit"); 
setAutoThreshold("Default dark");
//run("Threshold...");
setThreshold(1,255);
setOption("BlackBackground", true);
run("Convert to Mask");

selectWindow(poresfin);
saveAs("Tiff",poredir+origname+"_"+"Total_Pores.tif");
poresfin=getTitle();

selectWindow(pores);
close();

//Duplicate flattened pore image 

selectImage(poresfin); run("Duplicate...", "title=[CorTemp]");
cortemp="CorTemp";

selectImage(poresfin); run("Duplicate...", "title=[TrabTemp]");
trabtemp="TrabTemp";

//Analyze particles on flattened pore image

selectImage(poresfin);
run("Analyze Particles...", "display clear add");
run("Clear Results");

//Save new total pore ROIs

//roiManager("Deselect");
//roiManager("Save", wholedir+origname+"_"+"TotalPore_RoiSet.zip");

selectWindow(poresfin);
close();


//**Regional Subdivision****************************************************************************************************************

showStatus("!Regional Analysis...");

//User chooses rib or long bone 

Dialog.createNonBlocking("Cross-Section Type Selection");

Dialog.setInsets(5, 35, 0)
Dialog.addMessage("Select Cross-Section Type:", 14,"#7f0000");
Dialog.addRadioButtonGroup("", newArray("Rib","Long Bone"), 2, 1, "Rib");

Dialog.setInsets(0, 50, 0)
Dialog.addMessage("Draw Quadrants Using:");
Dialog.setInsets(0, 50, 0)
Dialog.addRadioButtonGroup("", newArray("Section Alignment With Image Borders","Section Major Axis"), 2, 1, "Section Alignment With Image Borders");

Dialog.show();

bonetype = Dialog.getRadioButton();
tilt = Dialog.getRadioButton();;

//Regional subdivision for long bones****************************************************************************************************

if (bonetype == "Long Bone"){

//Set scale for cross-sectional geometry according to user input in mm

run("Set Scale...", "distance=scale known=1 pixel=1 unit=mm global");

//Run BoneJ Slice Geometry 

selectImage(cortical);
run("Select None");
run("Remove Overlay");

selectImage(cortical);
run("Slice Geometry", "bone=unknown bone_min=1 bone_max=255 slope=0.0000 y_intercept=0");

selectImage(cortical);
getPixelSize(unit, pw, ph);

//Get row numbers of columns

selectWindow("Results");
text = getInfo();
lines = split(text, "\n");
columns = split(lines[0], "\t");

if (columns[0]==" ")
 columns[0]= "Number";

//Pull variables from BoneJ results table and divide by pixel size to obtain pixel coordinates

 cX= getResult(columns[5],0)/pw;
 cY = getResult(columns[6],0)/pw;

 th = abs(getResult(columns[10],0));
 rMin = getResult(columns[11],0)/pw;
 rMax = getResult(columns[12],0)/pw;
 thPi = th + PI / 2;

 //Loop to define tilt value for major axis quadrant subdivision 

if (tilt == "Section Major Axis"){

//Define major axis - this will be vertical for a long image, and horizontal for a short image

Majorx1 = floor(cX - cos(-th) * 1.5 * rMax);
Majory1 = floor(cY + sin(-th) * 1.5 * rMax);
Majorx2 = floor(cX + cos(-th) * 1.5 * rMax);
Majory2 = floor(cY - sin(-th) * 1.5 * rMax);

//Slope of line 

Major_m = (Majory1 - Majory2)/(Majorx1 - Majorx2);

//Because the coordinates are inverted (increase from top to bottom of frame), a  negative slope inclines to the right, and a positive slope inclines to the left
//Angle of line compared to vertical axis is tan angle = 1/m for a positive slope and -1/m for a negative slope
//Since a right incline is negative in slope and resulting angle, but rotation is clockwise, add a negative sign so that line rotates properly
//Multiply by 180/PI to convert from radians to degrees

Major_angle=-(atan(1/Major_m) * 180/PI);

//Define minor axis

Minorx1 = floor(cX - cos(thPi) * 2 * rMin);
Minory1 = floor(cY - sin(thPi) * 2 * rMin);
Minorx2 = floor(cX + cos(thPi) * 2 * rMin);
Minory2 = floor(cY + sin(thPi) * 2 * rMin);

//Slope of line 

Minor_m = (Minory1 - Minory2)/(Minorx1 - Minorx2);

//Because the coordinates are inverted (increase from top to bottom of frame), a  negative slope inclines to the right, and a positive slope inclines to the left
//Angle of line compared to vertical axis is tan angle = 1/m for a positive slope and -1/m for a negative slope
//Since a right incline is negative in slope and resulting angle, but rotation is clockwise, add a negative sign so that line rotates properly
//Multiply by 180/PI to convert from radians to degrees

Minor_angle=-(atan(1/Minor_m) * 180/PI);

//Detrmine coordinates for quadrant lines

//Get width and height of image 

w = getWidth();
h = getHeight();

//Length of the diagonal via pythagorean theorem 

d = sqrt((h*h) + (w*w));

//Difference between diagonal and height 

hdiff=((d-h)/2);

//Extend 5000 pixels further beyond diagonal

hboost=hdiff+5000;

//Clear slice geometry results as all values have been extracted 

run("Clear Results");

var x1 = "";
var y1 = "";
var x2 = "";
var y2 = "";
var Final_angle = "";

var linestop = "lineredo";

selectImage(cortical);
setBatchMode("show");
run("Select None");
run("Remove Overlay");

do {

//Erase ROIs 

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//For the adjusted centerline, make it the boosted diagonal length so it is sure to go out of frame by 1000 pixels even at the longest (diagonal) rotation

selectImage(cortical);
makeLine(cX, -hboost, cX, (h+hboost),30);

//Rotate line to major axis tilt angle 

run("Rotate...", "  angle=Major_angle");

Roi.setStrokeColor("red")
roiManager("Add");
roiManager("Show All without labels");

//Ask user to inspect regional subdivision

Dialog.createNonBlocking("Confirm Longest Axis");
Dialog.setInsets(5, 0, 5)
Dialog.addRadioButtonGroup("Are you satisfied with longest axis?",  newArray("Yes, Proceed", "No, Use Minor Axis"), 2, 1, "Yes, Proceed");

Dialog.show();

lineswap= Dialog.getRadioButton();

//If user chooses to proceed with major axis, exit do while loop

if (lineswap == "Yes, Proceed"){
	Final_angle = Major_angle;
	linestop = "linestop";
	}
else{linestop = "lineredo";}

//Swap names if chosen by user

if (lineswap == "No, Use Minor Axis"){

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//For the adjusted centerline, make it the boosted diagonal length so it is sure to go out of frame by 1000 pixels even at the longest (diagonal) rotation

selectImage(cortical);
makeLine(cX, -hboost, cX, (h+hboost),30);

//Rotate line to minor axis tilt angle 

run("Rotate...", "  angle=Minor_angle");

Roi.setStrokeColor("red")
roiManager("Add");
roiManager("Show All without labels");

//Ask user to inspect regional subdivision

Dialog.createNonBlocking("Confirm Longest Axis");
Dialog.setInsets(5, 0, 5)
Dialog.addRadioButtonGroup("Are you satisfied with longest axis?",  newArray("Yes, Proceed", "No, Use Major Axis"), 2, 1, "Yes, Proceed");

Dialog.show();

lineswapredo= Dialog.getRadioButton();

//If user chooses to proceed, exit do while loop
//If user chooses to revert to major axis, the loop will repeat

if (lineswapredo == "Yes, Proceed"){
	Final_angle = Minor_angle;
	linestop = "linestop";
	}
else{linestop = "lineredo";}

}

}


while (linestop=="lineredo");

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Hide total image again

selectImage(total);
setBatchMode("hide");

}

else{
	Final_angle = 0;
	run("Clear Results");
	};

//Then proceed from final angle - re-run extraction in case no tilt was selected

//Detrmine coordinates for quadrant lines

//Get width and height of image 

w = getWidth();
h = getHeight();

//Length of the diagonal via pythagorean theorem 

d = sqrt((h*h) + (w*w));

//Difference between diagonal and height 

hdiff=((d-h)/2);

//Extend 5000 pixels further beyond diagonal

hboost=hdiff+5000;

//For the adjusted centerline, make it the boosted diagonal length so it is sure to go out of frame by 1000 pixels even at the longest (diagonal) rotation

makeLine(cX, -hboost, cX, (h+hboost));

//Rotate line to tilt angle 

run("Rotate...", "  angle=Final_angle");


//Rotate major axis right 45 degrees to top of octant 2

run("Rotate...", "  angle=45");

//Get coordinates for opposing lines 2 and 4

getSelectionCoordinates(x, y);

o2x=x[0];
o2y=y[0];

o4x=x[1];
o4y=y[1];


//Rotate major axis right 90 degrees to top of octant 3

run("Rotate...", "  angle=90");

//Get coordinates for opposing lines 3 and 7

getSelectionCoordinates(x, y);

o3x=x[0];
o3y=y[0];

o1x=x[1];
o1y=y[1];

//Duplicate image for drawing all quadrants on the slice

selectImage(cortical); 
run("Select None");
run("Remove Overlay");
run("Duplicate...", "title=[Drawn]");

drawquad=getImageID();

//Change foreground to black to draw in black 

run("Colors...", "foreground=black background=black selection=red");
run("Line Width...", "line=1");

//Drawn octants on duplicate

selectImage(drawquad);
run("Select None");
run("Remove Overlay");

makePolygon(o1x,o1y,o2x,o2y,cX,cY);
run("Draw");
makePolygon(o3x,o3y,o4x,o4y,cX,cY);
run("Draw");
run("Select None");

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Add and combine in ROI manager 

makePolygon(o1x,o1y,o2x,o2y,cX,cY);
roiManager("Add");
makePolygon(o3x,o3y,o4x,o4y,cX,cY);
roiManager("Add");

roiManager("Select", newArray(0,1));
roiManager("XOR");
roiManager("Add");
roiManager("Delete");

run("From ROI Manager");

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}


//Change foreground back to white 

run("Colors...", "foreground=white background=black selection=red");

//Create folder subdivision for regional output 

regiondir=dir+"/Regions/";
File.makeDirectory(regiondir); 

selectImage(drawquad);
saveAs("TIFF", regiondir+origname+"_"+"DrawnQuadrants.tif");

drawquad=getImageID();
selectImage(drawquad);
close();

//Close the total area image

selectImage(total);
close();

//Clear outside Octant 1 on CA

selectImage(cortical); 
run("Select None");
run("Remove Overlay");

selectImage(cortical); 
run("Duplicate...", "title=[Quad1]");

quad1=getTitle();

selectImage(quad1); 

makePolygon(o1x,o1y,o2x,o2y,cX,cY);

setBackgroundColor(0, 0, 0);
run("Clear Outside");

setAutoThreshold("Default dark");
//run("Threshold...");
setThreshold(1,255);
setOption("BlackBackground", true);
run("Convert to Mask");

//Quadrant 1 Bounding
//Counts particle(s), if multiple ROIs combines them in a single ROI and deletes individual ROIs

run("Clear Results");

selectImage(quad1); 

run("Analyze Particles...", "display clear add");

var roicount = "";

roicount=roiManager("Count");

//Combine multiple ROIs for fragmented cortex region

if (roicount>1){

roiManager("show all without labels");

roiManager("Combine");

roiManager("Add");

newcount=roiManager("Count")-1;

deleteroi=Array.getSequence(newcount);

roiManager("Select", deleteroi);

roiManager("Delete");

roiManager("Select", 0);
}

//Or select single ROI for non-fragmented cortex region

else
{roiManager("Select", 0);}

run("Clear Results");

//Temporarily save ROI

roiManager("Select", 0);
roiManager("Save", regiondir+"TempQuad1.roi");

//Clear results 

run("Clear Results");

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Close quad 1 image 

selectImage(quad1);
close();

//Clear outside Quadrant 2 on CA

selectImage(cortical); 
run("Select None");
run("Remove Overlay");

selectImage(cortical); 
run("Duplicate...", "title=[Quad2]");

quad2=getTitle();

selectImage(quad2); 

makePolygon(o2x,o2y,o3x,o3y,cX,cY);

setBackgroundColor(0, 0, 0);
run("Clear Outside");

setAutoThreshold("Default dark");
//run("Threshold...");
setThreshold(1,255);
setOption("BlackBackground", true);
run("Convert to Mask");

//Quadrant 2 Bounding
//Counts particle(s), if multiple ROIs combines them in a single ROI and deletes individual ROIs

run("Clear Results");

selectImage(quad2); 

run("Analyze Particles...", "display clear add");

var roicount = "";

roicount=roiManager("Count");

//Combine multiple ROIs for fragmented cortex region

if (roicount>1){

roiManager("show all without labels");

roiManager("Combine");

roiManager("Add");

newcount=roiManager("Count")-1;

deleteroi=Array.getSequence(newcount);

roiManager("Select", deleteroi);

roiManager("Delete");

roiManager("Select", 0);
}

//Or select single ROI for non-fragmented cortex region

else
{roiManager("Select", 0);}

run("Clear Results");

//Temporarily save ROI

roiManager("Select", 0);
roiManager("Save", regiondir+"TempQuad2.roi");


//Clear results 

run("Clear Results");

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Close quad 2 image 

selectImage(quad2);
close();

//Clear outside Quadrant 3 on CA

selectImage(cortical); 
run("Select None");
run("Remove Overlay");

selectImage(cortical); 
run("Duplicate...", "title=[Quad3]");

quad3=getTitle();

selectImage(quad3); 

makePolygon(o3x,o3y,o4x,o4y,cX,cY);

setBackgroundColor(0, 0, 0);
run("Clear Outside");

setAutoThreshold("Default dark");
//run("Threshold...");
setThreshold(1,255);
setOption("BlackBackground", true);
run("Convert to Mask");

//Quadrant 3 Bounding
//Counts particle(s), if multiple ROIs combines them in a single ROI and deletes individual ROIs

run("Clear Results");

selectImage(quad3); 

run("Analyze Particles...", "display clear add");

var roicount = "";

roicount=roiManager("Count");

//Combine multiple ROIs for fragmented cortex region

if (roicount>1){

roiManager("show all without labels");

roiManager("Combine");

roiManager("Add");

newcount=roiManager("Count")-1;

deleteroi=Array.getSequence(newcount);

roiManager("Select", deleteroi);

roiManager("Delete");

roiManager("Select", 0);
}

//Or select single ROI for non-fragmented cortex region

else
{roiManager("Select", 0);}

run("Clear Results");

//Temporarily save ROI

roiManager("Select", 0);
roiManager("Save", regiondir+"TempQuad3.roi");


//Clear results 

run("Clear Results");

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}


//Close quad 3 image 

selectImage(quad3);
close();



//Clear outside Quadrant 4 on CA

selectImage(cortical); 
run("Select None");
run("Remove Overlay");

selectImage(cortical); 
run("Duplicate...", "title=[Quad4]");

quad4=getTitle();

selectImage(quad4); 

makePolygon(o4x,o4y,o1x,o1y,cX,cY);

setBackgroundColor(0, 0, 0);
run("Clear Outside");

setAutoThreshold("Default dark");
//run("Threshold...");
setThreshold(1,255);
setOption("BlackBackground", true);
run("Convert to Mask");

//Quadrant 4 Bounding
//Counts particle(s), if multiple ROIs combines them in a single ROI and deletes individual ROIs

run("Clear Results");

selectImage(quad4); 

run("Analyze Particles...", "display clear add");

var roicount = "";

roicount=roiManager("Count");

//Combine multiple ROIs for fragmented cortex region

if (roicount>1){

roiManager("show all without labels");

roiManager("Combine");

roiManager("Add");

newcount=roiManager("Count")-1;

deleteroi=Array.getSequence(newcount);

roiManager("Select", deleteroi);

roiManager("Delete");

roiManager("Select", 0);
}

//Or select single ROI for non-fragmented cortex region

else
{roiManager("Select", 0);}

run("Clear Results");

//Temporarily save ROI

roiManager("Select", 0);
roiManager("Save", regiondir+"TempQuad4.roi");


//Clear results 

run("Clear Results");

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Close quad 4 image 

selectImage(quad4);
close();

//Reopen ROIs

selectImage(cortical);
run("Select None");
run("Remove Overlay");

setBatchMode("show");

cortical = getTitle();

roiManager("Open", regiondir+"TempQuad1.roi"); 
roiManager("Open", regiondir+"TempQuad2.roi");
roiManager("Open", regiondir+"TempQuad3.roi"); 
roiManager("Open", regiondir+"TempQuad4.roi");  

//Use labels as names 

roiManager("UseNames", "true");

roiManager("Select", 0);
roiManager("Rename", "Quadrant_1");
roiManager("Set Fill Color", "red");

roiManager("Select", 1);
roiManager("Rename", "Quadrant_2");
roiManager("Set Fill Color", "blue");

roiManager("Select", 2);
roiManager("Rename", "Quadrant_3");
roiManager("Set Fill Color", "yellow");

roiManager("Select", 3);
roiManager("Rename", "Quadrant_4");
roiManager("Set Fill Color", "green");

selectImage(cortical);
roiManager("Show None");
roiManager("Show All with labels");


var regionstop = "regionredo";

//Do-while loop for user to rename quadrants

do {

//Dialog for quadrants set by anatomical orientation (directional display)

if (tilt == "Section Alignment With Image Borders"){

orientation=newArray("Anterior", "Posterior", "Medial", "Lateral", "Superior", "Inferior", "Custom");

Dialog.createNonBlocking("Quadrant Naming");
Dialog.addMessage("Choose anatomical orientations as they appear on the image");
Dialog.addMessage("");
Dialog.setInsets(5, 85, 5)
Dialog.addChoice("", orientation, "Anterior");
Dialog.setInsets(5, 0, 5)
Dialog.addChoice("Orientation:", orientation, "Lateral");
Dialog.addToSameRow();
Dialog.setInsets(5, 0, 5)
Dialog.addChoice("", orientation,"Medial");
Dialog.setInsets(5, 85, 5)
Dialog.addChoice("", orientation,"Posterior");
Dialog.show();

//Get user choices

regionname1 = Dialog.getChoice();
regionname4 = Dialog.getChoice();
regionname2 = Dialog.getChoice();
regionname3 = Dialog.getChoice();

}

//Dialog for quadrants set by major axis (no directional display)

else{

orientation=newArray("Anterior", "Posterior", "Medial", "Lateral", "Superior", "Inferior", "Custom");

Dialog.createNonBlocking("Quadrant Naming");
Dialog.addMessage("Choose anatomical orientations as they appear on the image");
Dialog.addMessage("");
Dialog.setInsets(5, 0, 5)
Dialog.addChoice("Quadrant 1 :", orientation, "Anterior");
Dialog.setInsets(5, 0, 5)
Dialog.addChoice("Quadrant 2 :", orientation, "Medial");
Dialog.setInsets(5, 0, 5)
Dialog.addChoice("Quadrant 3 :", orientation,"Posterior");
Dialog.setInsets(5, 0, 5)
Dialog.addChoice("Quadrant 4: ", orientation,"Lateral");
Dialog.show();

//Get user choices

regionname1 = Dialog.getChoice();
regionname2 = Dialog.getChoice();
regionname3 = Dialog.getChoice();
regionname4 = Dialog.getChoice();

	
}

//Loop for custom user-entered names

if(regionname1 == "Custom" || regionname2 == "Custom" || regionname3 == "Custom" || regionname4 == "Custom"){

loopcount=0;

Dialog.createNonBlocking("Custom Quadrant Name Entry");
Dialog.addMessage("Enter custom names for anatomical regions");

if(regionname1 == "Custom")
{Dialog.setInsets(0, 0, 0);
Dialog.addString("Quadrant 1 : ", "");
loopcount=loopcount+1;}
else
{Dialog.setInsets(0, 0, 0);
Dialog.addMessage("Quadrant 1 : " + regionname1);
}

if(regionname2 == "Custom")
{Dialog.setInsets(0, 0, 0);
Dialog.addString("Quadrant 2 :", "");
loopcount=loopcount+1;}
else
{Dialog.setInsets(0, 0, 0);
Dialog.addMessage("Quadrant 2 : " + regionname2);}

if(regionname3 == "Custom")
{Dialog.setInsets(0, 0, 0);
Dialog.addString("Quadrant 3 :", "");
loopcount=loopcount+1;}
else
{Dialog.setInsets(0, 0, 0);
Dialog.addMessage("Quadrant 3 : " +regionname3);}

if(regionname4 == "Custom")
{Dialog.setInsets(0, 0, 0);
Dialog.addString("Quadrant 4 :", "");
loopcount=loopcount+1;}
else
{Dialog.setInsets(0, 0, 0);
Dialog.addMessage("Quadrant 4 : " +regionname4);}

Dialog.show();

//Fill an array with the entered strings

loopout = newArray(loopcount);

for(i = 0; i<loopcount; i++){
looptemp = Dialog.getString();
loopout[i] = looptemp;
}

//Make an array of current names of regions

regionarray=newArray(regionname1,regionname2,regionname3,regionname4);

//Loop through the region name array and replace each next instance of "Custom" with the next instance of the custom string from loopout

loopoutpos=0;

for(i = 0; i<regionarray.length; i++){
	name=regionarray[i];
	if(name == "Custom"){
	regionarray[i] = loopout[loopoutpos];
	loopoutpos=loopoutpos+1;
	}
}

//Rename regions from the corrected regionarray

regionname1=regionarray[0];
regionname2=regionarray[1];
regionname3=regionarray[2];
regionname4=regionarray[3];

//Custom value correction end
}

//Rename ROIs

roiManager("Select", 0);
roiManager("Rename", regionname1);

roiManager("Select", 1);
roiManager("Rename", regionname2);

roiManager("Select", 2);
roiManager("Rename", regionname3);

roiManager("Select", 3);
roiManager("Rename", regionname4);

selectImage(cortical);
roiManager("Show None");
roiManager("Show All with labels");

//Ask user to inspect region names

Dialog.createNonBlocking("Confirm Quadrant Naming");
Dialog.setInsets(5, 0, 5)
Dialog.addRadioButtonGroup("Are you satisfied with quadrant labels?",  newArray("Yes, Proceed", "No, Rename Regions"), 2, 1, "Yes, Proceed");

Dialog.show();

regionredochoice= Dialog.getRadioButton();

if (regionredochoice == "No, Rename Regions"){

roiManager("Select", 0);
roiManager("Rename", "Quadrant_1");

roiManager("Select", 1);
roiManager("Rename", "Quadrant_2");

roiManager("Select", 2);
roiManager("Rename", "Quadrant_3");

roiManager("Select", 3);
roiManager("Rename", "Quadrant_4");

selectImage(cortical);
roiManager("Show None");
roiManager("Show All with labels");

regionredo = "regionredo";

}

else{

regionstop = "regionstop";

}
}

while (regionstop=="regionredo");


//Make subfolders for Regions folder for each region 

//Region 1 octant (Default: Superior) 

//regiondir1=regiondir+"/"+regionname1+"/";
//File.makeDirectory(regiondir1); 

//regiondir2=regiondir+"/"+regionname2+"/";
//File.makeDirectory(regiondir2); 

//regiondir3=regiondir+"/"+regionname3+"/";
//File.makeDirectory(regiondir3); 

//regiondir4=regiondir+"/"+regionname4+"/";
//File.makeDirectory(regiondir4); 

//Silence batch mode

setBatchMode("hide");

//Reset quadrant ROIs to blank color and fill and save to region subfolders, but do not close

roiManager("Select", 0);
roiManager("Set Color", "red");

roiManager("Select", 1);
roiManager("Set Color", "red");

roiManager("Select", 2);
roiManager("Set Color", "red");

roiManager("Select", 3);
roiManager("Set Color", "red");

selectImage(cortical);
roiManager("Show None");
roiManager("Show All with labels");

//Save ROIs to output folders

roiManager("Select",0);
roiManager("Save", regiondir+origname+"_"+regionname1+".roi");
roiManager("deselect");

roiManager("Select",1);
roiManager("Save", regiondir+origname+"_"+regionname2+".roi");
roiManager("deselect");

roiManager("Select",2);
roiManager("Save", regiondir+origname+"_"+regionname3+".roi");
roiManager("deselect");

roiManager("Select",3);
roiManager("Save", regiondir+origname+"_"+regionname4+".roi");
roiManager("deselect");


//Delete temp ROIs

File.delete(regiondir+"TempQuad1.roi");
File.delete(regiondir+"TempQuad2.roi");
File.delete(regiondir+"TempQuad3.roi");
File.delete(regiondir+"TempQuad4.roi");

selectWindow("Log");
run("Close");

//Clear results 

run("Clear Results");

//Set measurements to area only

run("Set Measurements...", "area redirect=None decimal=3");

//Measure each quadrant's CA

roiManager("Deselect");
roiManager("Measure");

var regionCA1 = "";
var regionCA2 = "";
var regionCA3 = "";
var regionCA4 = "";

var regionCA1 = getResult("Area",0);
var regionCA2 = getResult("Area",1);
var regionCA3 = getResult("Area",2);
var regionCA4 = getResult("Area",3);

//Close the cortical area image

selectImage(cortical);
close();

}

//Regional subdivision for ribs ******************************************************************************************************


if (bonetype == "Rib"){

//Set scale for cross-sectional geometry according to user input in um

run("Set Scale...", "distance=scale known=1 pixel=1 unit=mm global");

//Run BoneJ Slice Geometry 

selectImage(cortical);
run("Select None");
run("Remove Overlay");


selectImage(cortical);
run("Slice Geometry", "bone=unknown bone_min=1 bone_max=255 slope=0.0000 y_intercept=0");

selectImage(cortical);
getPixelSize(unit, pw, ph);

//Get image width and height - these are in pixels regardless of scale

w = getWidth();
h = getHeight();

//Get row numbers of columns

selectWindow("Results");
text = getInfo();
lines = split(text, "\n");
columns = split(lines[0], "\t");

if (columns[0]==" ")
 columns[0]= "Number";

//Pull variables from BoneJ results table and divide by pixel size
cX= getResult(columns[5],0)/pw;
cY = getResult(columns[6],0)/pw;

th = getResult(columns[10],0);
rMin = getResult(columns[11],0)/pw;
rMax = getResult(columns[12],0)/pw;
thPi = th + PI / 2;

//Define major axis

Majorx1 = floor(cX - cos(-th) * 1.5 * rMax);
Majory1 = floor(cY + sin(-th) * 1.5 * rMax);
Majorx2 = floor(cX + cos(-th) * 1.5 * rMax);
Majory2 = floor(cY - sin(-th) * 1.5 * rMax);

//drawLine(Majorx1, Majory1, Majorx2, Majory2);

//Define minor axis

Minorx1 = floor(cX - cos(thPi) * 2 * rMin);
Minory1 = floor(cY - sin(thPi) * 2 * rMin);
Minorx2 = floor(cX + cos(thPi) * 2 * rMin);
Minory2 = floor(cY + sin(thPi) * 2 * rMin);

var x1 = "";
var y1 = "";
var x2 = "";
var y2 = "";

var linestop = "lineredo";

selectImage(cortical);
setBatchMode("show");
run("Select None");
run("Remove Overlay");

do {

//Erase ROIs 

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}
	
//Draw the major axis 

selectImage(cortical);
makeLine(Majorx1,Majory1,Majorx2,Majory2,30);
Roi.setStrokeColor("red")
roiManager("Add");
roiManager("Show All without labels");

//Ask user to inspect regional subdivision

Dialog.createNonBlocking("Confirm Regional Subdivision");
Dialog.setInsets(5, 0, 5)
Dialog.addRadioButtonGroup("Are you satisfied with cutaneous / pleural subdivision?",  newArray("Yes, Proceed", "No, Use Minor Axis"), 2, 1, "Yes, Proceed");

Dialog.show();

lineswap= Dialog.getRadioButton();

//If user chooses to proceed with major axis, exit do while loop

if (lineswap == "Yes, Proceed"){
	x1 = Majorx1;
	y1= Majory1;
	x2 = Majorx2;
	y2 = Majory2;
	
	linestop = "linestop";
	}
else{linestop = "lineredo";}

//Swap names if chosen by user

if (lineswap == "No, Use Minor Axis"){

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Draw the minor axis 

selectImage(cortical);
makeLine(Minorx1,Minory1,Minorx2,Minory2,30);
Roi.setStrokeColor("red")
roiManager("Add");
roiManager("Show All without labels");

//Ask user to inspect regional subdivision

Dialog.createNonBlocking("Confirm Regional Subdivision");
Dialog.setInsets(5, 0, 5)
Dialog.addRadioButtonGroup("Are you satisfied with cutaneous / pleural subdivision?",  newArray("Yes, Proceed", "No, Use Major Axis"), 2, 1, "Yes, Proceed");

Dialog.show();

lineswapredo= Dialog.getRadioButton();

//If user chooses to proceed, exit do while loop
//If user chooses to revert to major axis, the loop will repeat

if (lineswapredo == "Yes, Proceed"){

	x1 = Minorx1;
	y1= Minory1;
	x2 = Minorx2;
	y2 = Minory2;
	linestop = "linestop";
	}
else{linestop = "lineredo";}

}

}


while (linestop=="lineredo");

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Hide cortical image again 

selectImage(cortical);
setBatchMode("hide");

//Duplicate image for drawing the major axis on the slice

selectImage(cortical);
run("Select None");
run("Remove Overlay");
run("Duplicate...", "title=[Drawn]");

drawhalf=getImageID();

//Change foreground to black

run("Colors...", "foreground=black background=black selection=red");
run("Line Width...", "line=1");

//Drawn major axis on the duplicate

selectImage(drawhalf);
run("Select None");
run("Remove Overlay");
makeLine(x1, y1, x2, y2);
run("Draw");
run("Add Selection...");
run("Select None");

//Change foreground back to white 

run("Colors...", "foreground=white background=black selection=red");

//Create folder subdivision for regional output 

regiondir=dir+"/Regions/";
File.makeDirectory(regiondir); 

selectImage(drawhalf);
saveAs("TIFF", regiondir+origname+"_"+"DrawnHalves.tif");

drawhalf=getImageID();
selectImage(drawhalf);
close();

//Close the total area image

selectImage(total);
close();


//Duplicate cortical image 

selectImage(cortical); run("Duplicate...", "title=[Region1]");
region1=getImageID();

selectImage(cortical); run("Duplicate...", "title=[Region2]");
region2=getImageID();

//Determine how the image is oriented with the long dimension of the rib 
//Vertically (height > width) or horizontally (width>height)
//Note: The axis drawn by the BoneJ macro is incorrect and extends the major axis beyond image bounds
//This corrected axis drawing will not match BoneJ macro axis output

setBackgroundColor(0, 0, 0);

setAutoThreshold("Default dark");
//run("Threshold...");
setThreshold(1,255);
setOption("BlackBackground", true);
run("Convert to Mask");


if (w>=h){

//Draw the top polygon on the horizontally oriented image
//x1,0 is the top left corner
//x1,0 is the top right corner

selectImage(region1);
makePolygon(x1,0,x2,0,x2,y2,x1,y1);
run("Clear Outside");
setAutoThreshold("Default dark");
//run("Threshold...");
setThreshold(1,255);
setOption("BlackBackground", true);
run("Convert to Mask");


//Draw the bottom polygon on the horizontally oriented image
//x1,h is the bottom left corner
//x2,h is the bottom right corner

selectImage(region2);
makePolygon(x1,h,x2,h,x2,y2,x1,y1);
run("Clear Outside");
setAutoThreshold("Default dark");
//run("Threshold...");
setThreshold(1,255);
setOption("BlackBackground", true);
run("Convert to Mask");


}
else 
{

//Draw the left polygon on the vertically oriented image
//0,y1 is upper left corner
//0,y2 is lower left corner

selectImage(region1);
makePolygon(0,y1,0,y2,x2,y2,x1,y1);
run("Clear Outside");
setAutoThreshold("Default dark");
//run("Threshold...");
setThreshold(1,255);
setOption("BlackBackground", true);
run("Convert to Mask");


//Draw the right polygon on the vertically oriented image
//w,y1 is upper right corner
//w,y2 is lower right corner

selectImage(region2);
makePolygon(w,y1,w,y2,x2,y2,x1,y1);
run("Clear Outside");
setAutoThreshold("Default dark");
//run("Threshold...");
setThreshold(1,255);
setOption("BlackBackground", true);
run("Convert to Mask");


}


//Clear results 

run("Clear Results");

//Clear ROI

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Region 1 measurement

//Analyze particles

selectImage(region1); 
run("Remove Overlay");
run("Select None");

run("Analyze Particles...", "display clear add");

//Combine multiple ROIs for fragmented cortex region

var roicount = "";

roicount=roiManager("Count");

if (roicount>1){

roiManager("show all without labels");

roiManager("Combine");

roiManager("Add");

newcount=roiManager("Count")-1;

deleteroi=Array.getSequence(newcount);

roiManager("Select", deleteroi);

roiManager("Delete");

roiManager("Select", 0);
}

//Or select single ROI for non-fragmented cortex region

else
{roiManager("Select", 0);}



//Temporarily save ROI

roiManager("Select", 0);
roiManager("Save", regiondir+"TempRegion1.roi");

//Clear results 

run("Clear Results");

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Close region 1 image 

selectImage(region1);
close();


//Region 2 measurement

//Analyze particles

selectImage(region2); 
run("Remove Overlay");
run("Select None");

run("Analyze Particles...", "display clear add");

//Combine multiple ROIs for fragmented cortex region

var roicount = "";

roicount=roiManager("Count");

if (roicount>1){

roiManager("show all without labels");

roiManager("Combine");

roiManager("Add");

newcount=roiManager("Count")-1;

deleteroi=Array.getSequence(newcount);

roiManager("Select", deleteroi);

roiManager("Delete");

roiManager("Select", 0);
}

//Or select single ROI for non-fragmented cortex region

else
{roiManager("Select", 0);}



//Temporarily save ROI

roiManager("Select", 0);
roiManager("Save", regiondir+"TempRegion2.roi");

//Clear results 

run("Clear Results");

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Close region 2 image 

selectImage(region2);
close();


//Reopen ROIs

selectImage(cortical);
run("Select None");
run("Remove Overlay");

setBatchMode("show");

cortical = getTitle();

roiManager("Open", regiondir+"TempRegion1.roi"); 
roiManager("Open", regiondir+"TempRegion2.roi");


////Guess cutaneous vs. pleural orientation 

//Set measurements to include area and shape descriptors

run("Set Measurements...", "area shape redirect=None decimal=3");

//Measure area and shape descriptors for each region

roiManager("Deselect");
roiManager("Measure");

var regionCA1 = "";
var regionCA2 = "";

var regionCA1 = getResult("Area",0);
Circ1 = getResult("Circ.",0);
//AR1 = getResult("AR",0);
//Round1 = getResult("Round",0);
//Solidity1 = getResult("Solidity",0);

var regionCA2 = getResult("Area",1);
Circ2 = getResult("Circ.",1);
//AR2 = getResult("AR",1);
//Round2 = getResult("Round",1);
//Solidity2 = getResult("Solidity",1);

//Test whether region 1 or region 2 has higher circularity (likely pleural)

var region1guess = "";
var region1alt = "";
var region1fin = "";
var region2guess = "";
var region2alt = "";
var region2fin = "";

if (Circ1>Circ2){
	region1guess = "Pleural";
	region1alt = "Cutaneous";
	region2guess = "Cutaneous";
	region2alt = "Pleural";
	}
else{
	region1guess = "Cutaneous";
	region1alt = "Pleural";
	region2guess = "Pleural";
	region2alt = "Cutaneous";
}


var regionstop = "regionredo";

do {

//Use labels as names 

roiManager("UseNames", "true");

//Rename ROIs based on guess

roiManager("Select", 0);
roiManager("Rename", region1guess);
roiManager("Set Fill Color", "red");

roiManager("Select", 1);
roiManager("Rename", region2guess);
roiManager("Set Fill Color", "blue");

selectImage(cortical);
roiManager("Show None");
roiManager("Show All with labels");

//Ask user to inspect region names

Dialog.createNonBlocking("Confirm Rib Region Naming");
Dialog.setInsets(5, 0, 5)
Dialog.addRadioButtonGroup("Are you satisfied with rib region labels?",  newArray("Yes, Proceed", "No, Swap Names"), 2, 1, "Yes, Proceed");

Dialog.show();

regionswap= Dialog.getRadioButton();

//If user chooses to proceed with guessed names, exit do while loop

if (regionswap == "Yes, Proceed"){
	region1fin = region1guess;
	region2fin = region2guess;
	regionstop = "regionstop";
	}
else{regionstop = "regionredo";}

//Swap names if chosen by user

if (regionswap == "No, Swap Names"){

roiManager("Select", 0);
roiManager("Rename", region1alt);

roiManager("Select", 1);
roiManager("Rename", region2alt);

selectImage(cortical);
roiManager("Show None");
roiManager("Show All with labels");

//Ask user to inspect region names

Dialog.createNonBlocking("Confirm Rib Region Naming");
Dialog.setInsets(5, 0, 5)
Dialog.addRadioButtonGroup("Are you satisfied with rib region labels?",  newArray("Yes, Proceed", "No, Swap Names"), 2, 1, "Yes, Proceed");

Dialog.show();

regionswapredo= Dialog.getRadioButton();

//If user chooses to proceed, exit do while loop
//If user chooses to revert the labels to the original, the loop will repeat

if (regionswapredo == "Yes, Proceed"){
	region1fin = region1alt;
	region2fin = region2alt;
	regionstop = "regionstop";}
else{regionstop = "regionredo";}

}

}


while (regionstop=="regionredo");

//Make subfolders for Regions folder for each region 

//Region 1 octant (Default: Superior) 

//regiondir1=regiondir+"/"+region1fin+"/";
//File.makeDirectory(regiondir1); 

//regiondir2=regiondir+"/"+region2fin+"/";
//File.makeDirectory(regiondir2); 

//Silence batch mode

setBatchMode("hide");

//Reset rib ROIs to blank color and fill and save to region subfolders, but do not close

roiManager("Select", 0);
roiManager("Set Color", "red");

roiManager("Select", 1);
roiManager("Set Color", "red");

selectImage(cortical);
roiManager("Show None");
roiManager("Show All with labels");

//Save ROIs to output folders

roiManager("Select",0);
roiManager("Save", regiondir+origname+"_"+region1fin+".roi");
roiManager("deselect");

roiManager("Select",1);
roiManager("Save", regiondir+origname+"_"+region2fin+".roi");
roiManager("deselect");

//Delete temp ROIs

File.delete(regiondir+"TempRegion1.roi");
File.delete(regiondir+"TempRegion2.roi");

selectWindow("Log");
run("Close");

//Clear results 

run("Clear Results");

//Close the cortical area image

selectImage(cortical);
close();

}

//Regional assignment******************************************************************************************
//This section opens each regional ROI and flattens it as a binary image
//The mean gray value is measured for each ROI on each region 
//ROIs fully within a region will be 255, and fully outside a region will be 0
//ROIs on the border between two regions will have the higher mean gray value in the region containing more of their area
//ROIs are labeled with their majority regional assignment and colorized in ROI Manager

function RegionMeasure(origpath,regionroipath,roipath,regionname) {

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Clear any past results 

run("Clear Results");

//Open the preprocessed DIC image as a template for the region mask

open(origpath); 

regionmask=getTitle();

//Clear origimg for region mask

selectImage(regionmask);
run("Select All");
setBackgroundColor(0, 0, 0);
run("Clear", "slice");

//Select and flatten region ROI 

roiManager("open", regionroipath);

selectImage(regionmask);
roiManager("Select", 0);
run("Fill", "slice");

selectImage(regionmask);
run("8-bit");
setAutoThreshold("Otsu");
//run("Threshold...");
setThreshold(1, 255);
setOption("BlackBackground", true);
run("Convert to Mask");

//Clear ROI manager

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Open border ROIs

roiManager("open", roipath);

//Renumber

for (i = 0; i<roiManager("count"); i++){
		roiManager("select", i);
		roiManager("rename", i)}

//Set measurements to mean gray values

run("Set Measurements...", "mean redirect=None decimal=3");

//Superimpose rois on regional mask 

selectImage(regionmask);

roiManager("Deselect");

roiManager("Measure");

//Print the mean gray value for each pore 

selectWindow("Results");

regionmean=Table.getColumn("Mean");

selectWindow("Regional Gray Mean");

Table.setColumn(regionname, regionmean);

selectImage(regionmask);
close();

}

//Create output table

Table.create("Regional Gray Mean");

regiongray="[Regional Gray Mean]"; 

selectWindow("Regional Gray Mean");
setLocation(screenWidth, screenHeight);

//Find mean gray value for each pore in each regional mask 

if (bonetype == "Long Bone"){

regionroipath1 = regiondir+origname+"_"+regionname1+".roi";
regionroipath2 = regiondir+origname+"_"+regionname2+".roi";
regionroipath3 = regiondir+origname+"_"+regionname3+".roi";
regionroipath4 = regiondir+origname+"_"+regionname4+".roi";


RegionMeasure(origpath,regionroipath1,roipath,regionname1);
RegionMeasure(origpath,regionroipath2,roipath,regionname2);
RegionMeasure(origpath,regionroipath3,roipath,regionname3);
RegionMeasure(origpath,regionroipath4,roipath,regionname4);

}

if (bonetype == "Rib"){

regionroipath1 = regiondir+origname+"_"+region1fin+".roi";
regionroipath2 = regiondir+origname+"_"+region2fin+".roi";

regionname1 = region1fin;
regionname2 = region2fin;

RegionMeasure(origpath,regionroipath1,roipath,regionname1);
RegionMeasure(origpath,regionroipath2,roipath,regionname2);

}

//Close DIC image

selectImage(origimg);
close();


//Clear any past results 

run("Clear Results");

//Add ROI row numbers

selectWindow("Regional Gray Mean");

roirow = Array.getSequence(roiManager("count"));

Table.setColumn("ROI", roirow);

//Assign region name based on max mean gray value in row 

roicount = roiManager("count");

regionnames = newArray(roicount);

for (i=0; i<roicount; i++){ 
 
selectWindow("Regional Gray Mean");

if (bonetype == "Long Bone"){

meanrow = newArray(Table.get(regionname1, i), Table.get(regionname2, i), Table.get(regionname3, i), Table.get(regionname4, i));

meanrownames = newArray(regionname1,regionname2,regionname3,regionname4);

//Last value [3] matches rank of max value
//Match this to the position of the row name 

meanrank = Array.rankPositions(meanrow);

maxmeanrank = meanrank[3];


}

if (bonetype == "Rib"){

meanrow = newArray(Table.get(regionname1, i), Table.get(regionname2, i));

meanrownames = newArray(regionname1,regionname2);

//Last value [1] matches rank of max value
//Match this to the position of the row name 

meanrank = Array.rankPositions(meanrow);

maxmeanrank = meanrank[1];

}

regionnames[i] = meanrownames[maxmeanrank];

}

Table.setColumn("Region", regionnames);


//Colorize ROIs based on region
//Also set ROI group based on region

roicount = roiManager("count");

if (bonetype == "Long Bone"){
	regioncolor = newArray("red", "yellow", "green", "blue");
	regiongroup = newArray(2, 7, 3, 1);
	}

if (bonetype == "Rib"){regioncolor = newArray("red", "blue");
	regiongroup = newArray(2, 1);
	}

for(i = 0; i<lengthOf(meanrownames); i++){

//Internal loop to change color if it matches a given region name

for(j = 0; j<roicount; j++){
	if (regionnames[j] == meanrownames[i]){
	roiManager("select", j);
	RoiManager.setGroup(0);
	RoiManager.setGroup(regiongroup[i]);
	roiManager("Set Fill Color", regioncolor[i]);
	};
}

}


//Save regionally colorized ROI set 

roiManager("Deselect");
roiManager("Save", poredir+origname+"_Total_Pores.zip");

//Save regional gray mean table 

selectWindow("Regional Gray Mean");
saveAs("Text", regiondir+origname+"_"+"Regional_Gray_Mean"+".csv");
run("Close");


//PORE MORPHOMETRIC ANALYSIS BY PORE TYPE ******************************************************************************************

showStatus("!Aggregate Pore Morphometric Analysis...");

//Make table for aggregate measurements 

summaryporemean="[Summary Pore Measurements]"; 
run("New... ", "name="+summaryporemean+" type=Table"); 
print(summaryporemean,"\\Headings:Image\tPore Type\tRegion\tCortical Area (um^2)\tCortical Area (mm^2)\t% Cortical Area\tBone Area (um^2)\tBone Area (mm^2)\t% Bone Area\tPercent Porosity (%)\tPore Density (1/um^2)\tPore Density (1/mm^2)\tTotal Pore Number\tTotal Pore Area (um^2)\tMean Pore Area (um^2)\tMean Pore Perimeter (um)\tMean Pore Ellipse Major Axis (um)\tMean Pore Ellipse Minor Axis (um)\tMean Pore Major Axis Angle (um)\tMean Pore Circularity\tMean Pore Max Feret Diameter (um)\tMean Pore Max Feret X (um)\tMean Pore Max Feret Y (um)\tMean Pore Feret Angle (degree)\tMean Min Pore Feret Diameter (um)\tMean Pore Aspect Ratio\tMean Pore Roundness\tMean Pore Solidity\t");

selectWindow("Summary Pore Measurements");
setLocation(screenWidth, screenHeight);

//TOTAL PORE MORPHOMETRY -----------------------------------------------------------------------------------------------------------------------------

//Convert CA and TA to um

CAum = CA*1000000;
TAum = TA*1000000;

if (bonetype == "Long Bone"){
regionCA1um = regionCA1*1000000;
regionCA2um = regionCA2*1000000;
regionCA3um = regionCA3*1000000;
regionCA4um = regionCA4*1000000;	

regionCAlist = newArray(regionCA1,regionCA2,regionCA3,regionCA4);
regionCAumlist = newArray(regionCA1um,regionCA2um,regionCA3um,regionCA4um);

}

if (bonetype == "Rib"){
regionCA1um = regionCA1*1000000;
regionCA2um = regionCA2*1000000;

regionCAlist = newArray(regionCA1,regionCA2);
regionCAumlist = newArray(regionCA1um,regionCA2um);

}


//Clear any current results 

run("Clear Results");

//Set scale in um according to user input

run("Set Scale...", "distance=scale_um known=1 pixel=1 unit=um global");

//Set measurements for pore systems 

run("Set Measurements...", "area min centroid perimeter fit shape feret's redirect=None decimal=3");

//Turn off row numbers for tables

setOption("ShowRowNumbers", false); 

//Measure all pore ROIs

roiManager("Deselect");
roiManager("Measure");

//Calculate number of rows before summary 

sumbegin = nResults;
Total_Number = sumbegin;


//Print individual pore output only if there is at least one total pore

if (Total_Number > 0){

//Make total pore output table

tot="[Individual Measurements]"; 
run("New... ", "name="+tot+" type=Table"); 
print(tot,"\\Headings:Pore\tRegion\tArea (um^2)\tCentroid X (um)\tCentroid Y (um)\tPerimeter (um)\tEllipse Major Axis (um)\tEllipse Minor Axis (um)\tMajor Axis Angle (degree)\tCircularity\tMax Feret Diameter (um)\tMax Feret X (um)\tMax Feret Y (um)\tFeret Angle (degree)\tMin Feret Diameter (um)\tAspect Ratio\tRoundness\tSolidity\t");

selectWindow("Results");

//Save values to total pore table

Group=Table.getColumn("Group");
Area=Table.getColumn("Area");
X=Table.getColumn("X");
Y=Table.getColumn("Y");
Perim=Table.getColumn("Perim.");
Major=Table.getColumn("Major");
Minor=Table.getColumn("Minor");
Angle=Table.getColumn("Angle");
Circ=Table.getColumn("Circ.");
Feret=Table.getColumn("Feret");
FeretX=Table.getColumn("FeretX");
FeretY=Table.getColumn("FeretY");
FeretAngle=Table.getColumn("FeretAngle");
MinFeret=Table.getColumn("MinFeret");
AR=Table.getColumn("AR");
Round=Table.getColumn("Round");
Solidity=Table.getColumn("Solidity");

//Change group designation 

poreregion = newArray(roicount);

if (bonetype == "Long Bone"){

for (j=0; j<Group.length; j++)
{
if (Group[j] == 2){poreregion[j] = regionname1;}
if (Group[j] == 7){poreregion[j] = regionname2;}
if (Group[j] == 3){poreregion[j] = regionname3;}
if (Group[j] == 1){poreregion[j] = regionname4;}
}
}


if (bonetype == "Rib"){
for (j=0; j<Group.length; j++)
{
if (Group[j] == 2){poreregion[j] = regionname1;}
if (Group[j] == 1){poreregion[j] = regionname2;}
}
}
	


for (j=0; j<nResults; j++)
      print(tot, j +"\t"+ poreregion[j] +"\t"+
      Area[j] +"\t"+ X[j] +"\t"+ Y[j] +"\t"+ 
	  Perim[j] +"\t"+ Major[j] +"\t"+ Minor[j] +"\t"+ Angle[j] +"\t"+ 
	  Circ[j]+"\t"+ Feret[j] +"\t"+ FeretX[j] +"\t"+ FeretY[j] +"\t"+ 
	  FeretAngle[j] +"\t"+ MinFeret[j] +"\t"+ AR[j] +"\t"+ Round[j] +"\t"+ Solidity[j] +"\t");

selectWindow("Individual Measurements");
saveAs("Text", poredir+origname+"_"+"Individual_Measurements"+".csv");
run("Close");

//Calculate percent porosity from TA before summarizing

//Sum pore areas

Total_Pore_Area = 0;

for (j=0; j<nResults; j++)
{Total_Pore_Area = Total_Pore_Area +  getResult("Area", j);}

Percent_Porosity= (Total_Pore_Area/CAum)*100;

//Calculate pore density per um from CA

Pore_Density_mm = Total_Number/CA;

Pore_Density = Total_Number/CAum;

//Calculate total porosity-corrected RCA for cross-sectional geometry table 

Total_Pore_Area_mm = Total_Pore_Area/1000000; 

CApore = CA-Total_Pore_Area_mm;

RCA_tot_pore=(CApore/TA)*100;

//Calculate percent cortical area

CA_percent = (CA/TA)*100;

//Calculate bone area 

BA = CA-Total_Pore_Area_mm;
BAum = CAum - Total_Pore_Area;
BA_percent = (BA/TA)*100;

//Summary for totals 

selectWindow("Results");

run("Summarize");

Mean_Area=Table.get("Area",sumbegin);
Mean_Perim=Table.get("Perim.",sumbegin);
Mean_Major=Table.get("Major",sumbegin);
Mean_Minor=Table.get("Minor",sumbegin);
Mean_Angle=Table.get("Angle",sumbegin);
Mean_Circ=Table.get("Circ.",sumbegin);
Mean_Feret=Table.get("Feret",sumbegin);
Mean_FeretX=Table.get("FeretX",sumbegin);
Mean_FeretY=Table.get("FeretY",sumbegin);
Mean_FeretAngle=Table.get("FeretAngle",sumbegin);
Mean_MinFeret=Table.get("MinFeret",sumbegin);
Mean_AR=Table.get("AR",sumbegin);
Mean_Round=Table.get("Round",sumbegin);
Mean_Solidity=Table.get("Solidity",sumbegin);


}


//If there are no total pores, print 0 for all values 

else 
{

Mean_Area=0;
Mean_Perim=0;
Mean_Major=0;
Mean_Minor=0;
Mean_Angle=0;
Mean_Circ=0;
Mean_Feret=0;
Mean_FeretX=0;
Mean_FeretY=0;
Mean_FeretAngle=0;
Mean_MinFeret=0;
Mean_AR=0;
Mean_Round=0;
Mean_Solidity=0;

Total_Pore_Area = 0;
Percent_Porosity= 0;
Pore_Density = 0;
Pore_Density_mm = 0;

//For no pores, bone area = cortical area values 

CA_percent = (CA/TA)*100;
BA = CA;
BAum = CAum;
BA_percent = CA_percent;


}


//Print to summary table 

regionname = "Total";

poretype="Total";

print(summaryporemean, origname +"\t"+ poretype +"\t"+ regionname +"\t"+ CAum +"\t"+ CA +"\t"+ CA_percent +"\t"+ BAum +"\t"+ BA +"\t"+ BA_percent +"\t"+ Percent_Porosity +"\t"+ Pore_Density +"\t"+ Pore_Density_mm +"\t"+ Total_Number +"\t"+
      Total_Pore_Area +"\t"+ Mean_Area +"\t"+ Mean_Perim +"\t"+ Mean_Major +"\t"+ Mean_Minor +"\t"+ Mean_Angle +"\t"+ 
	  Mean_Circ+"\t"+ Mean_Feret +"\t"+ Mean_FeretX +"\t"+ Mean_FeretY +"\t"+ 
	  Mean_FeretAngle +"\t"+ Mean_MinFeret +"\t"+ Mean_AR +"\t"+ Mean_Round +"\t"+ Mean_Solidity +"\t");


//Clear any current results 

run("Clear Results");

//Reopen total table 

open(poredir+origname+"_"+"Individual_Measurements.csv");

tottabname=File.name;

selectWindow(tottabname);
setLocation(screenWidth, screenHeight);

//Aggregate region lists

if (bonetype == "Long Bone"){regionlist = newArray(regionname1,regionname2,regionname3,regionname4);}

if (bonetype == "Rib"){regionlist = newArray(regionname1,regionname2);}
	
//Summary for regions - begin loop 

for (i=0; i<regionlist.length; i++){

currentregion = regionlist[i];

currentCA = regionCAlist[i];

currentCAum = regionCAumlist[i];

//Empty arrays from previous loop

var Mean_Area_array = newArray;
var Mean_Perim_array = newArray;
var Mean_Major_array = newArray;
var Mean_Minor_array = newArray;
var Mean_Angle_array = newArray;
var Mean_Circ_array = newArray;
var Mean_Feret_array = newArray;
var Mean_FeretX_array = newArray;
var Mean_FeretY_array = newArray;
var Mean_FeretAngle_array = newArray;
var Mean_MinFeret_array = newArray;
var Mean_AR_array = newArray;
var Mean_Round_array = newArray;
var Mean_Solidity_array = newArray;

//Start loop through table 

selectWindow(tottabname); 

totalcount = Table.size;

for (j=0; j<totalcount; j++){

selectWindow(tottabname);

tempregionname = Table.getString("Region", j);

//Get all values here

Temp_Mean_Area=Table.get("Area (um^2)",j);
Temp_Mean_Perim=Table.get("Perimeter (um)",j);
Temp_Mean_Major=Table.get("Ellipse Major Axis (um)",j);
Temp_Mean_Minor=Table.get("Ellipse Minor Axis (um)",j);
Temp_Mean_Angle=Table.get("Major Axis Angle (degree)",j);
Temp_Mean_Circ=Table.get("Circularity",j);
Temp_Mean_Feret=Table.get("Max Feret Diameter (um)",j);
Temp_Mean_FeretX=Table.get("Max Feret X (um)",j);
Temp_Mean_FeretY=Table.get("Max Feret Y (um)",j);
Temp_Mean_FeretAngle=Table.get("Feret Angle (degree)",j);
Temp_Mean_MinFeret=Table.get("Min Feret Diameter (um)",j);
Temp_Mean_AR=Table.get("Aspect Ratio",j);
Temp_Mean_Round=Table.get("Roundness",j);
Temp_Mean_Solidity=Table.get("Solidity",j);

//Assign these values to type arrays if they match region

if (tempregionname == currentregion){

Mean_Area_array = Array.concat(Mean_Area_array, Temp_Mean_Area);
Mean_Perim_array = Array.concat(Mean_Perim_array, Temp_Mean_Perim);
Mean_Major_array = Array.concat(Mean_Major_array, Temp_Mean_Major);
Mean_Minor_array = Array.concat(Mean_Minor_array, Temp_Mean_Minor);
Mean_Angle_array = Array.concat(Mean_Angle_array, Temp_Mean_Angle);
Mean_Circ_array = Array.concat(Mean_Circ_array, Temp_Mean_Circ);
Mean_Feret_array = Array.concat(Mean_Feret_array, Temp_Mean_Feret);
Mean_FeretX_array = Array.concat(Mean_FeretX_array, Temp_Mean_FeretX);
Mean_FeretY_array = Array.concat(Mean_FeretY_array, Temp_Mean_FeretY);
Mean_FeretAngle_array = Array.concat(Mean_FeretAngle_array, Temp_Mean_FeretAngle);
Mean_MinFeret_array = Array.concat(Mean_MinFeret_array, Temp_Mean_MinFeret);
Mean_AR_array = Array.concat(Mean_AR_array, Temp_Mean_AR);
Mean_Round_array = Array.concat(Mean_Round_array, Temp_Mean_Round);
Mean_Solidity_array = Array.concat(Mean_Solidity_array, Temp_Mean_Solidity);

}

//End loop for filling arrays

}


//Check whether any pores were assigned to the region 

Total_Number = Mean_Area_array.length; 

if (Total_Number > 0){

//Calculate means 

Array.getStatistics(Mean_Area_array, min, max, mean, std);
Mean_Area = (mean);

Array.getStatistics(Mean_Perim_array, min, max, mean, std);
Mean_Perim = (mean);

Array.getStatistics(Mean_Major_array, min, max, mean, std);
Mean_Major = (mean);

Array.getStatistics(Mean_Minor_array, min, max, mean, std);
Mean_Minor = (mean);

Array.getStatistics(Mean_Angle_array, min, max, mean, std);
Mean_Angle = (mean);

Array.getStatistics(Mean_Circ_array, min, max, mean, std);
Mean_Circ = (mean);

Array.getStatistics(Mean_Feret_array, min, max, mean, std);
Mean_Feret = (mean);

Array.getStatistics(Mean_FeretX_array, min, max, mean, std);
Mean_FeretX = (mean);

Array.getStatistics(Mean_FeretY_array, min, max, mean, std);
Mean_FeretY = (mean);

Array.getStatistics(Mean_FeretAngle_array, min, max, mean, std);
Mean_FeretAngle = (mean);

Array.getStatistics(Mean_MinFeret_array, min, max, mean, std);
Mean_MinFeret = (mean);

Array.getStatistics(Mean_AR_array, min, max, mean, std);
Mean_AR = (mean);

Array.getStatistics(Mean_Round_array, min, max, mean, std);
Mean_Round = (mean);

Array.getStatistics(Mean_Solidity_array, min, max, mean, std);
Mean_Solidity = (mean);

//Calculations 

//Sum pore areas

Total_Pore_Area = 0;

for (k=0; k<Mean_Area_array.length; k++)
{Total_Pore_Area = Total_Pore_Area + Mean_Area_array[k];}

Percent_Porosity= (Total_Pore_Area/currentCAum)*100;

//Calculate pore density per um from CA

Pore_Density_mm = Total_Number/currentCA;

Pore_Density = Total_Number/currentCAum;

//Calculate percent cortical area

currentCA_percent = (currentCA/TA)*100;

//Calculate bone area 

Total_Pore_Area_mm = Total_Pore_Area/1000000; 

BA = currentCA - Total_Pore_Area_mm;
BAum = currentCAum - Total_Pore_Area;
BA_percent = (BA/TA)*100;


}

//If there are no total pores in this region, print 0 for all values 

else 
{

Mean_Area=0;
Mean_Perim=0;
Mean_Major=0;
Mean_Minor=0;
Mean_Angle=0;
Mean_Circ=0;
Mean_Feret=0;
Mean_FeretX=0;
Mean_FeretY=0;
Mean_FeretAngle=0;
Mean_MinFeret=0;
Mean_AR=0;
Mean_Round=0;
Mean_Solidity=0;

Total_Pore_Area = 0;
Percent_Porosity= 0;
Pore_Density = 0;
Pore_Density_mm = 0;

//For no pores, bone area = cortical area values 

currentCA_percent = (currentCA/TA)*100;
BA = currentCA;
BAum = currentCAum;
BA_percent = currentCA_percent;


}

//Print to summary table 


print(summaryporemean, origname +"\t"+ poretype +"\t"+ currentregion +"\t"+ currentCAum +"\t"+ currentCA +"\t"+ currentCA_percent +"\t"+ BAum +"\t"+ BA +"\t"+ BA_percent +"\t"+ Percent_Porosity +"\t"+ Pore_Density +"\t"+ Pore_Density_mm +"\t"+ Total_Number +"\t"+
      Total_Pore_Area +"\t"+ Mean_Area +"\t"+ Mean_Perim +"\t"+ Mean_Major +"\t"+ Mean_Minor +"\t"+ Mean_Angle +"\t"+ 
	  Mean_Circ+"\t"+ Mean_Feret +"\t"+ Mean_FeretX +"\t"+ Mean_FeretY +"\t"+ 
	  Mean_FeretAngle +"\t"+ Mean_MinFeret +"\t"+ Mean_AR +"\t"+ Mean_Round +"\t"+ Mean_Solidity +"\t");


//End loop through regions
}


//PORE TYPE DIFFERENTIATION -------------------------------------------------------------------------------------------------------

showStatus("!Pore Type Differentiation...");

//Clear any past results

run("Clear Results");

//Set 32 bit for EDM map

run("Options...", "edm=32-bit");

//Remove scale to measure in pixels

selectImage(marrow);

run("Set Scale...", "distance=0 known=0 pixel=1 unit=pixel");

//Set measurements for pore systems 

run("Set Measurements...", "area min centroid perimeter fit shape feret's redirect=None decimal=3");

//Convert marrow to Euclidean Distance Map

selectImage(marrow); 
run("Distance Map"); 
//rename("EDM"); 
EDM=getTitle();

//Run measurements for the pore ROIs superimposed on the EDM of the marrow.
//Increasing minimum gray values on the EDM = increasing minimum distance from the marrow cavity
//Trabecularized if diameter (minferet) is greater than or equal to distance from the marrow (min gray value)

selectImage(EDM); 
roiManager("Show All");
roiManager("Measure");

//Close marrow and EDM image 

selectImage(marrow); 
close();

selectImage(EDM); 
close();

//Get min feret from results and convert to pixels

selectWindow("Results");

minferetpx=Table.getColumn("MinFeret");
mindistancepx=Table.getColumn("Min");

poretypearray = newArray(minferetpx.length);

for (i = 0; i < minferetpx.length; i++) {
	if(mindistancepx[i]<=minferetpx[i])
	{poretypearray[i] = "Trabecularized";}
	else{poretypearray[i] = "Cortical";}
}


//Add pore type to table 

selectWindow(tottabname);

Table.setColumn("Pore Type", poretypearray);

showStatus("!Cortical Pore Morphometric Analysis...");

//SUMMARY TABLE - CORTICAL PORES - TOTAL

//Clear any current results 

run("Clear Results");

//Empty arrays from previous loop

var Mean_Area_array = newArray;
var Mean_Perim_array = newArray;
var Mean_Major_array = newArray;
var Mean_Minor_array = newArray;
var Mean_Angle_array = newArray;
var Mean_Circ_array = newArray;
var Mean_Feret_array = newArray;
var Mean_FeretX_array = newArray;
var Mean_FeretY_array = newArray;
var Mean_FeretAngle_array = newArray;
var Mean_MinFeret_array = newArray;
var Mean_AR_array = newArray;
var Mean_Round_array = newArray;
var Mean_Solidity_array = newArray;

//Start loop through table 

selectWindow(tottabname); 

totalcount = Table.size;

for (i=0; i<totalcount; i++){

selectWindow(tottabname);

poretype= Table.getString("Pore Type", i);

//Get all values here

Temp_Mean_Area=Table.get("Area (um^2)",i);
Temp_Mean_Perim=Table.get("Perimeter (um)",i);
Temp_Mean_Major=Table.get("Ellipse Major Axis (um)",i);
Temp_Mean_Minor=Table.get("Ellipse Minor Axis (um)",i);
Temp_Mean_Angle=Table.get("Major Axis Angle (degree)",i);
Temp_Mean_Circ=Table.get("Circularity",i);
Temp_Mean_Feret=Table.get("Max Feret Diameter (um)",i);
Temp_Mean_FeretX=Table.get("Max Feret X (um)",i);
Temp_Mean_FeretY=Table.get("Max Feret Y (um)",i);
Temp_Mean_FeretAngle=Table.get("Feret Angle (degree)",i);
Temp_Mean_MinFeret=Table.get("Min Feret Diameter (um)",i);
Temp_Mean_AR=Table.get("Aspect Ratio",i);
Temp_Mean_Round=Table.get("Roundness",i);
Temp_Mean_Solidity=Table.get("Solidity",i);

//Assign these values to type arrays if they are cortical

if (poretype == "Cortical"){

Mean_Area_array = Array.concat(Mean_Area_array, Temp_Mean_Area);
Mean_Perim_array = Array.concat(Mean_Perim_array, Temp_Mean_Perim);
Mean_Major_array = Array.concat(Mean_Major_array, Temp_Mean_Major);
Mean_Minor_array = Array.concat(Mean_Minor_array, Temp_Mean_Minor);
Mean_Angle_array = Array.concat(Mean_Angle_array, Temp_Mean_Angle);
Mean_Circ_array = Array.concat(Mean_Circ_array, Temp_Mean_Circ);
Mean_Feret_array = Array.concat(Mean_Feret_array, Temp_Mean_Feret);
Mean_FeretX_array = Array.concat(Mean_FeretX_array, Temp_Mean_FeretX);
Mean_FeretY_array = Array.concat(Mean_FeretY_array, Temp_Mean_FeretY);
Mean_FeretAngle_array = Array.concat(Mean_FeretAngle_array, Temp_Mean_FeretAngle);
Mean_MinFeret_array = Array.concat(Mean_MinFeret_array, Temp_Mean_MinFeret);
Mean_AR_array = Array.concat(Mean_AR_array, Temp_Mean_AR);
Mean_Round_array = Array.concat(Mean_Round_array, Temp_Mean_Round);
Mean_Solidity_array = Array.concat(Mean_Solidity_array, Temp_Mean_Solidity);

}

//End loop for filling arrays

}


//Check whether any pores were assigned to be cortical

Total_Number = Mean_Area_array.length; 

if (Total_Number > 0){

//Calculate means 

Array.getStatistics(Mean_Area_array, min, max, mean, std);
Mean_Area = (mean);

Array.getStatistics(Mean_Perim_array, min, max, mean, std);
Mean_Perim = (mean);

Array.getStatistics(Mean_Major_array, min, max, mean, std);
Mean_Major = (mean);

Array.getStatistics(Mean_Minor_array, min, max, mean, std);
Mean_Minor = (mean);

Array.getStatistics(Mean_Angle_array, min, max, mean, std);
Mean_Angle = (mean);

Array.getStatistics(Mean_Circ_array, min, max, mean, std);
Mean_Circ = (mean);

Array.getStatistics(Mean_Feret_array, min, max, mean, std);
Mean_Feret = (mean);

Array.getStatistics(Mean_FeretX_array, min, max, mean, std);
Mean_FeretX = (mean);

Array.getStatistics(Mean_FeretY_array, min, max, mean, std);
Mean_FeretY = (mean);

Array.getStatistics(Mean_FeretAngle_array, min, max, mean, std);
Mean_FeretAngle = (mean);

Array.getStatistics(Mean_MinFeret_array, min, max, mean, std);
Mean_MinFeret = (mean);

Array.getStatistics(Mean_AR_array, min, max, mean, std);
Mean_AR = (mean);

Array.getStatistics(Mean_Round_array, min, max, mean, std);
Mean_Round = (mean);

Array.getStatistics(Mean_Solidity_array, min, max, mean, std);
Mean_Solidity = (mean);

//Calculations 

//Sum pore areas

Total_Pore_Area = 0;

for (k=0; k<Mean_Area_array.length; k++)
{Total_Pore_Area = Total_Pore_Area + Mean_Area_array[k];}

Percent_Porosity= (Total_Pore_Area/CAum)*100;

//Calculate pore density per um from CA

Pore_Density_mm = Total_Number/CA;

Pore_Density = Total_Number/CAum;

//Calculate cortical porosity-corrected RCA for cross-sectional geometry table 

Total_Pore_Area_mm = Total_Pore_Area/1000000; 

CApore = CA-Total_Pore_Area_mm;

RCA_cor_pore=(CApore/TA)*100;

//Calculate percent cortical area

CA_percent = (CA/TA)*100;

//Calculate bone area 

BA = CA-Total_Pore_Area_mm;
BAum = CAum - Total_Pore_Area;
BA_percent = (BA/TA)*100;

}

//If there are no cortical pores, print 0 for all values 

else 
{

Mean_Area=0;
Mean_Perim=0;
Mean_Major=0;
Mean_Minor=0;
Mean_Angle=0;
Mean_Circ=0;
Mean_Feret=0;
Mean_FeretX=0;
Mean_FeretY=0;
Mean_FeretAngle=0;
Mean_MinFeret=0;
Mean_AR=0;
Mean_Round=0;
Mean_Solidity=0;

Total_Pore_Area = 0;
Percent_Porosity= 0;
Pore_Density = 0;
Pore_Density_mm = 0;

RCA_cor_pore=0;

//For no pores, bone area = cortical area values 

CA_percent = (CA/TA)*100;
BA = CA;
BAum = CAum;
BA_percent = CA_percent;

}

//Print to summary table 

poretype = "Cortical";

currentregion = "Total";

print(summaryporemean, origname +"\t"+ poretype +"\t"+ regionname +"\t"+ CAum +"\t"+ CA +"\t"+ CA_percent +"\t"+ BAum +"\t"+ BA +"\t"+ BA_percent +"\t"+ Percent_Porosity +"\t"+ Pore_Density +"\t"+ Pore_Density_mm +"\t"+ Total_Number +"\t"+
      Total_Pore_Area +"\t"+ Mean_Area +"\t"+ Mean_Perim +"\t"+ Mean_Major +"\t"+ Mean_Minor +"\t"+ Mean_Angle +"\t"+ 
	  Mean_Circ+"\t"+ Mean_Feret +"\t"+ Mean_FeretX +"\t"+ Mean_FeretY +"\t"+ 
	  Mean_FeretAngle +"\t"+ Mean_MinFeret +"\t"+ Mean_AR +"\t"+ Mean_Round +"\t"+ Mean_Solidity +"\t");



showStatus("!Regional Cortical Pore Morphometric Analysis...");

//SUMMARY TABLE - CORTICAL PORES - REGIONAL 

//Clear any current results 

run("Clear Results");

	
//Summary for regions - begin loop 

for (i=0; i<regionlist.length; i++){

currentregion = regionlist[i];

currentCA = regionCAlist[i];

currentCAum = regionCAumlist[i];

//Empty arrays from previous loop

var Mean_Area_array = newArray;
var Mean_Perim_array = newArray;
var Mean_Major_array = newArray;
var Mean_Minor_array = newArray;
var Mean_Angle_array = newArray;
var Mean_Circ_array = newArray;
var Mean_Feret_array = newArray;
var Mean_FeretX_array = newArray;
var Mean_FeretY_array = newArray;
var Mean_FeretAngle_array = newArray;
var Mean_MinFeret_array = newArray;
var Mean_AR_array = newArray;
var Mean_Round_array = newArray;
var Mean_Solidity_array = newArray;

//Start loop through table 

selectWindow(tottabname); 

totalcount = Table.size;

for (j=0; j<totalcount; j++){

selectWindow(tottabname);

tempregionname = Table.getString("Region", j);
poretype= Table.getString("Pore Type", j);

//Get all values here

Temp_Mean_Area=Table.get("Area (um^2)",j);
Temp_Mean_Perim=Table.get("Perimeter (um)",j);
Temp_Mean_Major=Table.get("Ellipse Major Axis (um)",j);
Temp_Mean_Minor=Table.get("Ellipse Minor Axis (um)",j);
Temp_Mean_Angle=Table.get("Major Axis Angle (degree)",j);
Temp_Mean_Circ=Table.get("Circularity",j);
Temp_Mean_Feret=Table.get("Max Feret Diameter (um)",j);
Temp_Mean_FeretX=Table.get("Max Feret X (um)",j);
Temp_Mean_FeretY=Table.get("Max Feret Y (um)",j);
Temp_Mean_FeretAngle=Table.get("Feret Angle (degree)",j);
Temp_Mean_MinFeret=Table.get("Min Feret Diameter (um)",j);
Temp_Mean_AR=Table.get("Aspect Ratio",j);
Temp_Mean_Round=Table.get("Roundness",j);
Temp_Mean_Solidity=Table.get("Solidity",j);

//Assign these values to type arrays if they match region and are cortical

if (tempregionname == currentregion && poretype == "Cortical"){

Mean_Area_array = Array.concat(Mean_Area_array, Temp_Mean_Area);
Mean_Perim_array = Array.concat(Mean_Perim_array, Temp_Mean_Perim);
Mean_Major_array = Array.concat(Mean_Major_array, Temp_Mean_Major);
Mean_Minor_array = Array.concat(Mean_Minor_array, Temp_Mean_Minor);
Mean_Angle_array = Array.concat(Mean_Angle_array, Temp_Mean_Angle);
Mean_Circ_array = Array.concat(Mean_Circ_array, Temp_Mean_Circ);
Mean_Feret_array = Array.concat(Mean_Feret_array, Temp_Mean_Feret);
Mean_FeretX_array = Array.concat(Mean_FeretX_array, Temp_Mean_FeretX);
Mean_FeretY_array = Array.concat(Mean_FeretY_array, Temp_Mean_FeretY);
Mean_FeretAngle_array = Array.concat(Mean_FeretAngle_array, Temp_Mean_FeretAngle);
Mean_MinFeret_array = Array.concat(Mean_MinFeret_array, Temp_Mean_MinFeret);
Mean_AR_array = Array.concat(Mean_AR_array, Temp_Mean_AR);
Mean_Round_array = Array.concat(Mean_Round_array, Temp_Mean_Round);
Mean_Solidity_array = Array.concat(Mean_Solidity_array, Temp_Mean_Solidity);

}

//End loop for filling arrays

}


//Check whether any pores were assigned to the region 

Total_Number = Mean_Area_array.length; 

if (Total_Number > 0){

//Calculate means 

Array.getStatistics(Mean_Area_array, min, max, mean, std);
Mean_Area = (mean);

Array.getStatistics(Mean_Perim_array, min, max, mean, std);
Mean_Perim = (mean);

Array.getStatistics(Mean_Major_array, min, max, mean, std);
Mean_Major = (mean);

Array.getStatistics(Mean_Minor_array, min, max, mean, std);
Mean_Minor = (mean);

Array.getStatistics(Mean_Angle_array, min, max, mean, std);
Mean_Angle = (mean);

Array.getStatistics(Mean_Circ_array, min, max, mean, std);
Mean_Circ = (mean);

Array.getStatistics(Mean_Feret_array, min, max, mean, std);
Mean_Feret = (mean);

Array.getStatistics(Mean_FeretX_array, min, max, mean, std);
Mean_FeretX = (mean);

Array.getStatistics(Mean_FeretY_array, min, max, mean, std);
Mean_FeretY = (mean);

Array.getStatistics(Mean_FeretAngle_array, min, max, mean, std);
Mean_FeretAngle = (mean);

Array.getStatistics(Mean_MinFeret_array, min, max, mean, std);
Mean_MinFeret = (mean);

Array.getStatistics(Mean_AR_array, min, max, mean, std);
Mean_AR = (mean);

Array.getStatistics(Mean_Round_array, min, max, mean, std);
Mean_Round = (mean);

Array.getStatistics(Mean_Solidity_array, min, max, mean, std);
Mean_Solidity = (mean);

//Calculations 

//Sum pore areas

Total_Pore_Area = 0;

for (k=0; k<Mean_Area_array.length; k++)
{Total_Pore_Area = Total_Pore_Area + Mean_Area_array[k];}

Percent_Porosity= (Total_Pore_Area/currentCAum)*100;

//Calculate pore density per um from CA

Pore_Density_mm = Total_Number/currentCA;

Pore_Density = Total_Number/currentCAum;

//Calculate percent cortical area

currentCA_percent = (currentCA/TA)*100;

//Calculate bone area 

Total_Pore_Area_mm = Total_Pore_Area/1000000; 

BA = currentCA - Total_Pore_Area_mm;
BAum = currentCAum - Total_Pore_Area;
BA_percent = (BA/TA)*100;

}

//If there are no total pores in this region, print 0 for all values 

else 
{

Mean_Area=0;
Mean_Perim=0;
Mean_Major=0;
Mean_Minor=0;
Mean_Angle=0;
Mean_Circ=0;
Mean_Feret=0;
Mean_FeretX=0;
Mean_FeretY=0;
Mean_FeretAngle=0;
Mean_MinFeret=0;
Mean_AR=0;
Mean_Round=0;
Mean_Solidity=0;

Total_Pore_Area = 0;
Percent_Porosity= 0;
Pore_Density = 0;
Pore_Density_mm = 0;

//For no pores, bone area = cortical area values 

currentCA_percent = (currentCA/TA)*100;
BA = currentCA;
BAum = currentCAum;
BA_percent = currentCA_percent;

}

//Print to summary table 

poretype = "Cortical";

print(summaryporemean, origname +"\t"+ poretype +"\t"+ currentregion +"\t"+ currentCAum +"\t"+ currentCA +"\t"+ currentCA_percent +"\t"+ BAum +"\t"+ BA +"\t"+ BA_percent +"\t"+ Percent_Porosity +"\t"+ Pore_Density +"\t"+ Pore_Density_mm +"\t"+ Total_Number +"\t"+
      Total_Pore_Area +"\t"+ Mean_Area +"\t"+ Mean_Perim +"\t"+ Mean_Major +"\t"+ Mean_Minor +"\t"+ Mean_Angle +"\t"+ 
	  Mean_Circ+"\t"+ Mean_Feret +"\t"+ Mean_FeretX +"\t"+ Mean_FeretY +"\t"+ 
	  Mean_FeretAngle +"\t"+ Mean_MinFeret +"\t"+ Mean_AR +"\t"+ Mean_Round +"\t"+ Mean_Solidity +"\t");



//End loop through regions
}

showStatus("!Trabecularized Pore Morphometric Analysis...");

//SUMMARY TABLE - TRABECULARIZED PORES - TOTAL

//Clear any current results 

run("Clear Results");

//Empty arrays from previous loop

var Mean_Area_array = newArray;
var Mean_Perim_array = newArray;
var Mean_Major_array = newArray;
var Mean_Minor_array = newArray;
var Mean_Angle_array = newArray;
var Mean_Circ_array = newArray;
var Mean_Feret_array = newArray;
var Mean_FeretX_array = newArray;
var Mean_FeretY_array = newArray;
var Mean_FeretAngle_array = newArray;
var Mean_MinFeret_array = newArray;
var Mean_AR_array = newArray;
var Mean_Round_array = newArray;
var Mean_Solidity_array = newArray;

//Start loop through table 

selectWindow(tottabname); 

totalcount = Table.size;

for (i=0; i<totalcount; i++){

selectWindow(tottabname);

poretype= Table.getString("Pore Type", i);

//Get all values here

Temp_Mean_Area=Table.get("Area (um^2)",i);
Temp_Mean_Perim=Table.get("Perimeter (um)",i);
Temp_Mean_Major=Table.get("Ellipse Major Axis (um)",i);
Temp_Mean_Minor=Table.get("Ellipse Minor Axis (um)",i);
Temp_Mean_Angle=Table.get("Major Axis Angle (degree)",i);
Temp_Mean_Circ=Table.get("Circularity",i);
Temp_Mean_Feret=Table.get("Max Feret Diameter (um)",i);
Temp_Mean_FeretX=Table.get("Max Feret X (um)",i);
Temp_Mean_FeretY=Table.get("Max Feret Y (um)",i);
Temp_Mean_FeretAngle=Table.get("Feret Angle (degree)",i);
Temp_Mean_MinFeret=Table.get("Min Feret Diameter (um)",i);
Temp_Mean_AR=Table.get("Aspect Ratio",i);
Temp_Mean_Round=Table.get("Roundness",i);
Temp_Mean_Solidity=Table.get("Solidity",i);

//Assign these values to type arrays if they are cortical

if (poretype == "Trabecularized"){

Mean_Area_array = Array.concat(Mean_Area_array, Temp_Mean_Area);
Mean_Perim_array = Array.concat(Mean_Perim_array, Temp_Mean_Perim);
Mean_Major_array = Array.concat(Mean_Major_array, Temp_Mean_Major);
Mean_Minor_array = Array.concat(Mean_Minor_array, Temp_Mean_Minor);
Mean_Angle_array = Array.concat(Mean_Angle_array, Temp_Mean_Angle);
Mean_Circ_array = Array.concat(Mean_Circ_array, Temp_Mean_Circ);
Mean_Feret_array = Array.concat(Mean_Feret_array, Temp_Mean_Feret);
Mean_FeretX_array = Array.concat(Mean_FeretX_array, Temp_Mean_FeretX);
Mean_FeretY_array = Array.concat(Mean_FeretY_array, Temp_Mean_FeretY);
Mean_FeretAngle_array = Array.concat(Mean_FeretAngle_array, Temp_Mean_FeretAngle);
Mean_MinFeret_array = Array.concat(Mean_MinFeret_array, Temp_Mean_MinFeret);
Mean_AR_array = Array.concat(Mean_AR_array, Temp_Mean_AR);
Mean_Round_array = Array.concat(Mean_Round_array, Temp_Mean_Round);
Mean_Solidity_array = Array.concat(Mean_Solidity_array, Temp_Mean_Solidity);

}

//End loop for filling arrays

}


//Check whether any pores were assigned to be cortical

Total_Number = Mean_Area_array.length; 

if (Total_Number > 0){

//Calculate means 

Array.getStatistics(Mean_Area_array, min, max, mean, std);
Mean_Area = (mean);

Array.getStatistics(Mean_Perim_array, min, max, mean, std);
Mean_Perim = (mean);

Array.getStatistics(Mean_Major_array, min, max, mean, std);
Mean_Major = (mean);

Array.getStatistics(Mean_Minor_array, min, max, mean, std);
Mean_Minor = (mean);

Array.getStatistics(Mean_Angle_array, min, max, mean, std);
Mean_Angle = (mean);

Array.getStatistics(Mean_Circ_array, min, max, mean, std);
Mean_Circ = (mean);

Array.getStatistics(Mean_Feret_array, min, max, mean, std);
Mean_Feret = (mean);

Array.getStatistics(Mean_FeretX_array, min, max, mean, std);
Mean_FeretX = (mean);

Array.getStatistics(Mean_FeretY_array, min, max, mean, std);
Mean_FeretY = (mean);

Array.getStatistics(Mean_FeretAngle_array, min, max, mean, std);
Mean_FeretAngle = (mean);

Array.getStatistics(Mean_MinFeret_array, min, max, mean, std);
Mean_MinFeret = (mean);

Array.getStatistics(Mean_AR_array, min, max, mean, std);
Mean_AR = (mean);

Array.getStatistics(Mean_Round_array, min, max, mean, std);
Mean_Round = (mean);

Array.getStatistics(Mean_Solidity_array, min, max, mean, std);
Mean_Solidity = (mean);

//Calculations 

//Sum pore areas

Total_Pore_Area = 0;

for (k=0; k<Mean_Area_array.length; k++)
{Total_Pore_Area = Total_Pore_Area + Mean_Area_array[k];}

Percent_Porosity= (Total_Pore_Area/CAum)*100;

//Calculate pore density per um from CA

Pore_Density_mm = Total_Number/CA;

Pore_Density = Total_Number/CAum;

//Calculate trabecularized porosity-corrected RCA for cross-sectional geometry table 

Total_Pore_Area_mm = Total_Pore_Area/1000000; 

CApore = CA-Total_Pore_Area_mm;

RCA_trab_pore=(CApore/TA)*100;

//Calculate percent cortical area

CA_percent = (CA/TA)*100;

//Calculate bone area 

BA = CA-Total_Pore_Area_mm;
BAum = CAum - Total_Pore_Area;
BA_percent = (BA/TA)*100;

}

//If there are no total pores in this region, print 0 for all values 

else 
{

Mean_Area=0;
Mean_Perim=0;
Mean_Major=0;
Mean_Minor=0;
Mean_Angle=0;
Mean_Circ=0;
Mean_Feret=0;
Mean_FeretX=0;
Mean_FeretY=0;
Mean_FeretAngle=0;
Mean_MinFeret=0;
Mean_AR=0;
Mean_Round=0;
Mean_Solidity=0;

Total_Pore_Area = 0;
Percent_Porosity= 0;
Pore_Density = 0;
Pore_Density_mm = 0;

RCA_trab_pore = 0;


//For no pores, bone area = cortical area values 

CA_percent = (CA/TA)*100;
BA = CA;
BAum = CAum;
BA_percent = CA_percent;



}

//Print to summary table 

poretype = "Trabecularized";

currentregion = "Total";

print(summaryporemean, origname +"\t"+ poretype +"\t"+ regionname +"\t"+ CAum +"\t"+ CA +"\t"+ CA_percent +"\t"+ BAum +"\t"+ BA +"\t"+ BA_percent +"\t"+ Percent_Porosity +"\t"+ Pore_Density +"\t"+ Pore_Density_mm +"\t"+ Total_Number +"\t"+
      Total_Pore_Area +"\t"+ Mean_Area +"\t"+ Mean_Perim +"\t"+ Mean_Major +"\t"+ Mean_Minor +"\t"+ Mean_Angle +"\t"+ 
	  Mean_Circ+"\t"+ Mean_Feret +"\t"+ Mean_FeretX +"\t"+ Mean_FeretY +"\t"+ 
	  Mean_FeretAngle +"\t"+ Mean_MinFeret +"\t"+ Mean_AR +"\t"+ Mean_Round +"\t"+ Mean_Solidity +"\t");



//SUMMARY TABLE - TRABECULARIZED PORES - REGIONAL 

showStatus("!Regional Trabecularized Pore Morphometric Analysis...");

//Clear any current results 

run("Clear Results");

	
//Summary for regions - begin loop 

for (i=0; i<regionlist.length; i++){

currentregion = regionlist[i];

currentCA = regionCAlist[i];

currentCAum = regionCAumlist[i];

//Empty arrays from previous loop

var Mean_Area_array = newArray;
var Mean_Perim_array = newArray;
var Mean_Major_array = newArray;
var Mean_Minor_array = newArray;
var Mean_Angle_array = newArray;
var Mean_Circ_array = newArray;
var Mean_Feret_array = newArray;
var Mean_FeretX_array = newArray;
var Mean_FeretY_array = newArray;
var Mean_FeretAngle_array = newArray;
var Mean_MinFeret_array = newArray;
var Mean_AR_array = newArray;
var Mean_Round_array = newArray;
var Mean_Solidity_array = newArray;

//Start loop through table 

selectWindow(tottabname); 

totalcount = Table.size;

for (j=0; j<totalcount; j++){

selectWindow(tottabname);

tempregionname = Table.getString("Region", j);
poretype= Table.getString("Pore Type", j);

//Get all values here

Temp_Mean_Area=Table.get("Area (um^2)",j);
Temp_Mean_Perim=Table.get("Perimeter (um)",j);
Temp_Mean_Major=Table.get("Ellipse Major Axis (um)",j);
Temp_Mean_Minor=Table.get("Ellipse Minor Axis (um)",j);
Temp_Mean_Angle=Table.get("Major Axis Angle (degree)",j);
Temp_Mean_Circ=Table.get("Circularity",j);
Temp_Mean_Feret=Table.get("Max Feret Diameter (um)",j);
Temp_Mean_FeretX=Table.get("Max Feret X (um)",j);
Temp_Mean_FeretY=Table.get("Max Feret Y (um)",j);
Temp_Mean_FeretAngle=Table.get("Feret Angle (degree)",j);
Temp_Mean_MinFeret=Table.get("Min Feret Diameter (um)",j);
Temp_Mean_AR=Table.get("Aspect Ratio",j);
Temp_Mean_Round=Table.get("Roundness",j);
Temp_Mean_Solidity=Table.get("Solidity",j);

//Assign these values to type arrays if they match region and are cortical

if (tempregionname == currentregion && poretype == "Trabecularized"){

Mean_Area_array = Array.concat(Mean_Area_array, Temp_Mean_Area);
Mean_Perim_array = Array.concat(Mean_Perim_array, Temp_Mean_Perim);
Mean_Major_array = Array.concat(Mean_Major_array, Temp_Mean_Major);
Mean_Minor_array = Array.concat(Mean_Minor_array, Temp_Mean_Minor);
Mean_Angle_array = Array.concat(Mean_Angle_array, Temp_Mean_Angle);
Mean_Circ_array = Array.concat(Mean_Circ_array, Temp_Mean_Circ);
Mean_Feret_array = Array.concat(Mean_Feret_array, Temp_Mean_Feret);
Mean_FeretX_array = Array.concat(Mean_FeretX_array, Temp_Mean_FeretX);
Mean_FeretY_array = Array.concat(Mean_FeretY_array, Temp_Mean_FeretY);
Mean_FeretAngle_array = Array.concat(Mean_FeretAngle_array, Temp_Mean_FeretAngle);
Mean_MinFeret_array = Array.concat(Mean_MinFeret_array, Temp_Mean_MinFeret);
Mean_AR_array = Array.concat(Mean_AR_array, Temp_Mean_AR);
Mean_Round_array = Array.concat(Mean_Round_array, Temp_Mean_Round);
Mean_Solidity_array = Array.concat(Mean_Solidity_array, Temp_Mean_Solidity);

}

//End loop for filling arrays

}


//Check whether any pores were assigned to the region 

Total_Number = Mean_Area_array.length; 

if (Total_Number > 0){

//Calculate means 

Array.getStatistics(Mean_Area_array, min, max, mean, std);
Mean_Area = (mean);

Array.getStatistics(Mean_Perim_array, min, max, mean, std);
Mean_Perim = (mean);

Array.getStatistics(Mean_Major_array, min, max, mean, std);
Mean_Major = (mean);

Array.getStatistics(Mean_Minor_array, min, max, mean, std);
Mean_Minor = (mean);

Array.getStatistics(Mean_Angle_array, min, max, mean, std);
Mean_Angle = (mean);

Array.getStatistics(Mean_Circ_array, min, max, mean, std);
Mean_Circ = (mean);

Array.getStatistics(Mean_Feret_array, min, max, mean, std);
Mean_Feret = (mean);

Array.getStatistics(Mean_FeretX_array, min, max, mean, std);
Mean_FeretX = (mean);

Array.getStatistics(Mean_FeretY_array, min, max, mean, std);
Mean_FeretY = (mean);

Array.getStatistics(Mean_FeretAngle_array, min, max, mean, std);
Mean_FeretAngle = (mean);

Array.getStatistics(Mean_MinFeret_array, min, max, mean, std);
Mean_MinFeret = (mean);

Array.getStatistics(Mean_AR_array, min, max, mean, std);
Mean_AR = (mean);

Array.getStatistics(Mean_Round_array, min, max, mean, std);
Mean_Round = (mean);

Array.getStatistics(Mean_Solidity_array, min, max, mean, std);
Mean_Solidity = (mean);

//Calculations 

//Sum pore areas

Total_Pore_Area = 0;

for (k=0; k<Mean_Area_array.length; k++)
{Total_Pore_Area = Total_Pore_Area + Mean_Area_array[k];}

Percent_Porosity= (Total_Pore_Area/currentCAum)*100;

//Calculate pore density per um from CA

Pore_Density_mm = Total_Number/currentCA;

Pore_Density = Total_Number/currentCAum;


//Calculate percent cortical area

currentCA_percent = (currentCA/TA)*100;

//Calculate bone area 

Total_Pore_Area_mm = Total_Pore_Area/1000000; 

BA = currentCA - Total_Pore_Area_mm;
BAum = currentCAum - Total_Pore_Area;
BA_percent = (BA/TA)*100;


}

//If there are no total pores in this region, print 0 for all values 

else 
{

Mean_Area=0;
Mean_Perim=0;
Mean_Major=0;
Mean_Minor=0;
Mean_Angle=0;
Mean_Circ=0;
Mean_Feret=0;
Mean_FeretX=0;
Mean_FeretY=0;
Mean_FeretAngle=0;
Mean_MinFeret=0;
Mean_AR=0;
Mean_Round=0;
Mean_Solidity=0;

Total_Pore_Area = 0;
Percent_Porosity= 0;
Pore_Density = 0;
Pore_Density_mm = 0;

//For no pores, bone area = cortical area values 

currentCA_percent = (currentCA/TA)*100;
BA = currentCA;
BAum = currentCAum;
BA_percent = currentCA_percent;

}

//Print to summary table 

poretype = "Trabecularized";

print(summaryporemean, origname +"\t"+ poretype +"\t"+ currentregion +"\t"+ currentCAum +"\t"+ currentCA +"\t"+ currentCA_percent +"\t"+ BAum +"\t"+ BA +"\t"+ BA_percent +"\t"+ Percent_Porosity +"\t"+ Pore_Density +"\t"+ Pore_Density_mm +"\t"+ Total_Number +"\t"+
      Total_Pore_Area +"\t"+ Mean_Area +"\t"+ Mean_Perim +"\t"+ Mean_Major +"\t"+ Mean_Minor +"\t"+ Mean_Angle +"\t"+ 
	  Mean_Circ+"\t"+ Mean_Feret +"\t"+ Mean_FeretX +"\t"+ Mean_FeretY +"\t"+ 
	  Mean_FeretAngle +"\t"+ Mean_MinFeret +"\t"+ Mean_AR +"\t"+ Mean_Round +"\t"+ Mean_Solidity +"\t");

//End loop through regions
}

showStatus("!Saving Tables...");

//Print to cross-sectional geometry table 

print(rca,origname+"\t"+TA+"\t"+MA+"\t"+CA+"\t"+RCA+"\t"+RCA_tot_pore+"\t"+RCA_cor_pore+"\t"+RCA_trab_pore+"\t"+PerMA+"\t"+Para+"\t"+Imin+"\t"+Imax+"\t"+Zpol+"\t");

//Save output table

selectWindow("Relative Cortical Area");
saveAs("Text", geomdir+origname+"_"+"Cross_Sectional_Geometry"+".csv");

//Save individual pore measurements table 

selectWindow(tottabname); 
saveAs("Text", poredir+origname+"_"+"Individual_Measurements"+".csv");
//run("Close");

//Save summary table 

selectWindow("Summary Pore Measurements");
saveAs("Text", poredir+origname+"_"+"Summary_Measurements"+".csv");
run("Close");

//Save ROI sets and TIFF images for cortical and trabecularized pores---------------------------------------------------------------

showStatus("!Saving ROI Sets and Images...");

selectWindow(tottabname);

poretypeout = Table.getColumn("Pore Type");

poretypeindexcor = newArray;
poretypeindextrab = newArray;

for (i = 0; i < poretypeout.length; i++) {
	if (poretypeout[i] == "Cortical"){poretypeindexcor  = Array.concat(poretypeindexcor, i);}
	else {poretypeindextrab  = Array.concat(poretypeindextrab, i);}
}


//Array.print(poretypeindexcor);
//Array.print(poretypeindextrab);

//Save cortical pores

if (poretypeindexcor.length > 0){

roiManager("deselect");
roiManager("Select", poretypeindexcor);
roiManager("Save Selected", poredir+origname+"_Cortical_Pores.zip");
}

//Save trabecularized pores

if (poretypeindextrab.length > 0){

roiManager("deselect");
roiManager("Select", poretypeindextrab);
roiManager("Save Selected", poredir+origname+"_Trabecularized_Pores.zip");

}

//Clear ROI manager 

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//If both cortical and trabecularized pores have pores, then save images individiually 

if (poretypeindexcor.length > 0 && poretypeindextrab.length > 0){

//Flatten black cortical pores to get white image of trabecilarized pores

corpath = poredir+origname+"_Cortical_Pores.zip";
roiManager("open", corpath);

selectWindow(trabtemp);
run("Colors...", "foreground=black background=black selection=red");
roiManager("deselect");
roiManager("Fill");
roiManager("Show None");
run("Flatten");
selectWindow("TrabTemp-1");
rename("Trabecularized Pores"); 
trabfin="Trabecularized Pores";
run("Colors...", "foreground=white background=black selection=red");

selectWindow(trabfin);
run("8-bit"); 
setAutoThreshold("Default dark");
setThreshold(1,255);
setOption("BlackBackground", true);
run("Convert to Mask");

saveAs("Tiff",poredir+origname+"_Trabecularized_Pores.tif");
trabfin=getTitle();

//Clear ROI manager 

finalcount=roiManager("count")-1;

if(finalcount>=0){
roiManager("Deselect");
roiManager("Delete");
}

//Flatten black trabecularized pores to get white image of cortical pores

trabpath = poredir+origname+"_Trabecularized_Pores.zip";
roiManager("open", trabpath);

selectWindow(cortemp);
run("Colors...", "foreground=black background=black selection=red");
roiManager("deselect");
roiManager("Fill");
roiManager("Show None");
run("Flatten");
selectWindow("CorTemp-1");
rename("Cortical Pores"); 
corfin="Cortical Pores";
run("Colors...", "foreground=white background=black selection=red");

selectWindow(corfin);
run("8-bit"); 
setAutoThreshold("Default dark");
setThreshold(1,255);
setOption("BlackBackground", true);
run("Convert to Mask");

saveAs("Tiff",poredir+origname+"_Cortical_Pores.tif");
corfin=getTitle();

}

//Do not save any images if there are no cortical pores or trabecularized pores individually 


//Close any remaining images

while (nImages>0) { 
          selectImage(nImages); 
          close(); 
}

//Close any remaining windows

list = getList("window.titles");
     for (i=0; i<list.length; i++){
     winame = list[i];
      selectWindow(winame);
     run("Close");
     }

//Clear any past results

run("Clear Results");

showStatus("!Pore Analysis Complete!");

}



//-----------------------------------------------------------------------------------------------------------------------------------------

function About() {

	if(isOpen("Log")==1) { selectWindow("Log"); run("Close"); }
	Dialog.create("About");
	Dialog.addMessage("Copyright (C) 2022 Mary E. Cole / PoreExtractor 2D \n \n Redistribution and use in source and binary forms, with or without \n modification, are permitted provided that the following conditions are met:\n \n 1. Redistributions of source code must retain the above copyright notice, \n this list of conditions and the following disclaimer.\n \n 2. Redistributions in binary form must reproduce the above copyright \nnotice, this list of conditions and the following disclaimer in the \ndocumentation and/or other materials provided with the distribution. \n \n 3. Neither the name of PoreExtractor 2D nor the names of its \n contributors may be used to endorse or promote products derived from \nthis software without specific prior written permission.\n \n This program is free software: you can redistribute it and/or modify \n it under the terms of the GNU General Public License as published by \nthe Free Software Foundation, either version 3 of the License, or \n any later version. \n \n This program is distributed in the hope that it will be useful, \n but WITHOUT ANY WARRANTY; without even the implied warranty of \n MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the \n GNU General Public License for more details.");
	Dialog.addHelp("http://www.gnu.org/licenses/");
	Dialog.show();
	
}


function Cite() {

	if(isOpen("Log")==1) { selectWindow("Log"); run("Close"); }
	 Dialog.createNonBlocking("Cite");
	 Dialog.addMessage("Mary E. Cole, Samuel D. Stout, Amanda M. Agnew. 2021.\nPore Extractor 2D: An ImageJ plugin for identification, classification, and regional characterization of cortical pores on histological bone images.\nAmerican Journal of Physical Anthropology 174(S71): 20-21.");
	 Dialog.addMessage("\nSelect OK to copy citation to system clipboard.");
	 Dialog.show();  
	 String.copy("Mary E. Cole, Samuel D. Stout, Amanda M. Agnew. 2021. Pore Extractor 2D: An ImageJ plugin for identification, classification, and regional characterization of cortical pores on histological bone images. American Journal of Physical Anthropology 174(S71): 20-21.");
}