/**
  * discrete_histogram_entropy_.ijm
  * 
  * Calculates the discrete histogram entropy 
  * for the result of the Directionality plugin
  *   
  * (c) 2019, INSERM
  * 
  * written by Volker Baecker at Montpellier Ressources Imagerie (www.mri.cnrs.fr)
  *
  * USAGE:
  *  Run the Directionality plugin (see https://imagej.net/Directionality) on your image
  *  Select "Display table" in the dialog. Activate the results table containing for each direction
  *  the normalized frequency of the data and the fit. Run the macro. The discrete histogram entropy 
  *  is written to the log window.
  *  
  *  The macro is available on git-hub:
  *  https://github.com/MontpellierRessourcesImagerie/imagej_macros_and_scripts/blob/master/volker/macros/discrete_histogram_entropy_.ijm
  *  
  */

_NBINS = 256;
_DELTA = 0.001;
_VALUE_COLUMN = "bin start";
_COUNT_COLUMN = "count";
_AUTO = true;

entropy = autoHistogramEntropy();
exit();

macro "Discrete Histogram Entropy (f5)" {
	autoHistogramEntropy();
}

macro "Discrete Histogram Entropy (f5) Action Tool - CeefD10CfccL2030CccfL7080CfdeDc0CfccDd0CfdeDe0D01Ce56D11Cb46L2131Ce9aD41C44eD61C43cL7181C77dD91CfccDb1Ce56Lc1e1CfccDf1D02Cb46D12Ce11L2232Ce56D42CfdeD52C43cD62C44eL7282C43cD92CccfDa2Ce9aDb2Ce11Lc2e2Ce9aDf2CfccD03Cb46D13Ce11L2333Ce9aD43CeefD53Cb46D63C11dL7383C43cD93CeefDa3Ce9aDb3Cb46Dc3Ce11Dd3Cb46De3Ce9aDf3CfccD14Ce56D24Ce9aD34CfdeD44CeefD64C44eD74C77dD84CccfD94Ce9aDc4Ce56Dd4Ce9aDe4CfdeL2535CeefD75CccfD85CfdeDc5CfccDd5CfdeDe5C77dD16C43cL2636CccfD46D66C43cL7686C77dD96CfccDb6Ce56Dc6Cb46Dd6Ce56De6CfccDf6CccfD07C11dD17C44eD27C11dD37C77dD47CeefD57C43cD67C44eL7787C43cD97CccfDa7Ce9aDb7Ce11Lc7e7Ce9aDf7CccfD08C43cD18C11dL2838C77dD48D68C11dL7888C43cD98CfdeDa8Ce9aDb8Ce56Dc8Ce11Dd8Cb46De8Ce9aDf8CccfD19Cb46D29C77dD39CeefD49CccfD69C77dD79Cb46D89CccfD99CfdeDb9Ce9aDc9Ce56Dd9Ce9aDe9CccfD2aCeefD3aCfdeL7a8aCeefDcaCfdeDdaCeefDeaC44eD1bC43cD2bCb46D3bCccfD4bCfccD6bCe56D7bCb46D8bCe9aD9bCeefDbbC77dDcbC43cDdbC77dDebCccfD0cC11dD1cC44eD2cC11dD3cC77dD4cCfdeD5cCe56D6cCe11L7c8cCb46D9cCfccDacC44eDbcC11dDccC44eDdcC11dDecC44eDfcCccfD0dC11dD1dC44eD2dC11dD3dC77dD4dCeefD5dCe56D6dCe11L7d8dCb46D9dCfccDadC44eDbdC11dDcdC44eDddC11dDedC44eDfdD1eC43cD2eCb46D3eCccfD4eCfccD6eCe56L7e8eCe9aD9eC77dDceC43cDdeC77dDeeCeefL2f3fD8f" {
	autoHistogramEntropy();
}

function autoHistogramEntropy() {
	winType = getInfo("window.type");
	title = getInfo("window.title");
	if (winType=='Image') {
		entropy = imageEntropy();
		print("Entropy of " + title + ": "+entropy);
		return entropy;
	} 
	if (indexOf(title, 'Directionality histograms')>=0) {
		entropy = directonalityEntropy();
		return entropy;
	}
	valueColumn = _VALUE_COLUMN;
	countColumn = _COUNT_COLUMN;
	if (_AUTO) {
		countColumn = 'count';
		headings = Table.headings;
		valuesColumn = "index";
		if (indexOf(headings, 'bin start')>=0) {
			valuesColumn = "bin start";
		}
	}
	entropy = histogramEntropy(title, valuesColumn, countColumn);
	print("Entropy of " + title + ": "+entropy);
	return entropy;
}

function histogramEntropy(title, columnValues, columnCounts) {
	values = Table.getColumn(columnValues, title);
	counts = Table.getColumn(columnCounts, title);
	sum = sumArray(counts);
	diffFromOne = abs(sum-1);
	if (diffFromOne>_DELTA) {
		normalizeArray(counts);
	}
	entropy = discreteHistogramEntropy(values, counts);	
	return entropy;
}

function imageEntropy() {
	imageTitle = getTitle();
	getHistogram(values, counts, _NBINS);
	sum = sumArray(counts);
	diffFromOne = abs(sum-1);
	if (diffFromOne>_DELTA) {
		normalizeArray(counts);
	}
	entropy = discreteHistogramEntropy(values, counts);
	return entropy;
}

function normalizeArray(anArray) {
	sum = sumArray(anArray);
	for (i = 0; i < anArray.length; i++) {
		anArray[i] /= sum;
	}
}

function sumArray(anArray) {
	sum = 0;
	for (i = 0; i < anArray.length; i++) {
		sum += anArray[i];
	}
	return sum;
}

function directonalityEntropy() {
	content = getInfo("window.contents");
	lines = split(content, "\n");
	line0 = split(lines[0],"\t");
	imageTitle = line0[1];
	directions = newArray(lines.length-1);
	frequencies = newArray(lines.length-1);
	for (i = 1; i < lines.length; i++) {
		line =  split(lines[i],"\t");
		directions[i-1] = parseFloat(line[0]);
		frequencies[i-1] = parseFloat(line[1]);
	}	
	entropy = discreteHistogramEntropy(directions, frequencies);
	print("Entropy of " + imageTitle + ": "+entropy);
	return entropy;
}

function discreteHistogramEntropy(bins, frequencies) {

	binWidth = abs(bins[0]-bins[1]);
	
	entropy = 0;
	for(i=0; i<bins.length; i++) {
		if (frequencies[i]>0)
			entropy += frequencies[i]*log(frequencies[i]/binWidth);
	}
	
	entropy *= -1;
	return entropy;
}

function entropyFilter(radius) {
	imgWidth = getWidth();
	imgHeight = getHeight();
	imageID = getImageID();
	run("Duplicate...", " ");
	run("32-bit");
	outID = getImageID();
	width = (2 * radius) + 1;
	setBatchMode(true);
	for (i = 0; i < imgWidth; i++) {
		for (j = 0; j < imgHeight; j++) {
			selectImage(imageID);
			makeRectangle(i-radius, j-radius, width, width);
			entropy = imageEntropy();
			selectImage(outID);
			setPixel(i, j, entropy);
		}
	}
	setBatchMode(false);
	selectImage(imageID);
	run("Select None");
	selectImage(outID);
}