#pragma rtGlobals=1		// Use modern global access method.

// ================================================================
//
// 	Multilevel Thresholding based on a "maximized between-class variance" criterion
// 	originally due to N. Otsu (IEEE Trans System Man Cybernetics, 1979)
//
// 	This implementation follows the article "A Fast Algorithm for Multilevel Thresholding"
// 	by P.S. Liao, T.S. Chen, and T.C. Chung (=LCC), J Inf Sci Eng 17, 713-727 (2001)
//
// 	Implementation for IgorPro 6.1+  by  W. Harneit (Berlin, Germany)
//
// 	2010/01/27	-	version 0.99
//
// ================================================================
//
//	usage: 	1. Display a grayscale image with NewImage
//			2. Bring it to the front
//			3. On the command line, execute 	• MultiThresh("", n, 0)
//			    where n is the desired number of threshold levels
//
//	hints:	n=1 and n=2 are pretty fast, 		< 1 s
//			n=3 may take some time, and 		~ 7.5 s
//			n=4 is really slow.				~ 750 s
//		Times determined on a (single processor) 2.67 GHz / 1.0 GB WinXP machine,
//		with kNumBins=256 for image "lena" (taken from "Image Processing Tutorial.pxp",
//		to be found in "Igor Pro Folder:Learning Aids:Tutorials:") 
//		(use Browse Experiment button in Data Browser, available in Menu "Data")
//
//		Setting kNumBins < kMaxGrayLevels speeds up things considerably,
//		but it also leads to (interesting?) artefacts (try kNumBins = 64 on "MRI")
//
//	to do list:	• put all temporary data into a package data folder
//				• investigate artefacts when kNumBins < kMaxGrayLevels
//				• nicer user interface (panel with global parameters)
//				• better behaviour when called with specified image name
//
// ================================================================

constant kMaxThresholdLevels = 4, kMaxGrayLevels = 256, kNumBins = 256

// main entry point
function MultiThresh(imageName, nThresholds, forceNewTable)
string imageName
variable nThresholds, forceNewTable

	if( strlen(imageName) == 0 )
		imageName = StringFromList(0, ImageNameList("",";"))
	endif

	wave imageWave = ImageNameToWaveRef("", imageName)
	if( !WaveExists(imageWave) )
		abort "image not found"
	endif
	
	string tableName = imageName+"_LCCtable"
	if( !WaveExists($tableName) || forceNewTable )
		print "precomputing LCC tables..."
		PrecomputeTables(imageWave, tableName)
	endif
	
	print "finding optimal thresholds..."
	FindThresholds(tableName, nThresholds)
	
	print "creating transformed image..."
	ApplyThresholds(imageWave, imageName+"_Thresh"+num2str(nThresholds))
end

// create new image with multi-level threshold
// thresholds are expected in ThresholdWave = {0, k1, ..., kn, kMaxGrayLevel}
function ApplyThresholds(imageWave, newImageName)
wave imageWave
string newImageName
	Wave tWave = ThresholdWave		// created by FindThresholds

	Make/B/U/O/N=(kMaxGrayLevels) LookupTable	// create a lookup table
	variable n, nGrayLevels = numpnts(tWave)-1
	for( n = 0; n < nGrayLevels; n += 1 )
		LookupTable[tWave[n]+1,tWave[n+1]] = (n+1) / nGrayLevels * (kMaxGrayLevels-1)
	endfor
	
	Duplicate/O imageWave, $newImageName
	Wave nWave = $newImageName
	nWave = LookupTable[nWave]		// apply lookup transform
	
	NewImage/K=1 nWave			// display result in new image
	AutoPositionWindow $WinName(0,1)
end

// precompute some tables according to LCC
function PrecomputeTables(imageWave, tableName)
wave imageWave
string tableName
	Make/D/O/N=(kNumBins) LCC0, LCC1
	Make/D/O/N=(kNumBins-1,kNumBins-1) ptab, stab, htab
	Histogram/B={-0.5,1,kNumBins} imageWave, LCC0	// centered histogram of gray levels in image
	LCC0[0]=0						// zap zero values
	LCC1 = LCC0 * p					// first moment
	Integrate LCC0					// integrated (or cumulative) histogram is ...
	ptab[0][] = LCC0[q+1]				// ... first row (here: column) in LCC's P table
	ptab[1,kNumBins-1][] = ptab[0][q] - ptab[0][p-1]	// the other rows (columns)
	Integrate LCC1					// integrated first moment is ...
	stab[0][] = LCC1[q+1]				// ... first row (here: column) in LCC's S table
	stab[1,kNumBins-1][] = stab[0][q] - stab[0][p-1]	// the other rows (columns)
	htab = stab/ptab*stab				// S and P table yield LCC's H table
	Duplicate/O htab $tableName		// store for later re-use
	killwaves LCC0, LCC1, stab, ptab, htab
end

// this function just dispatches to the functions below (+ pre/post-processing)
function FindThresholds(tableName, nThresholds)
string tableName
variable nThresholds
	variable nLevels = limit(floor(nThresholds), 1, kMaxThresholdLevels)
	if( nLevels != nThresholds )
		print "number of threshold levels set to", nLevels
	endif
	Make/O/N=(nLevels+2) ThresholdWave = 0
	Execute "FindThresh"+num2str(nLevels)+"("+tableName+")"
	print ThresholdWave
	ThresholdWave[nLevels+1] = kNumBins-1
end

function FindThresh1(Tab)
wave Tab
	wave ThresholdWave
	variable k1, sdev, maxsdev = 0
	for( k1 = 1; k1 < kNumBins-1; k1 += 1 )
		sdev = Tab[1][k1] + Tab[k1+1][kNumBins-1]
		if( sdev > maxsdev )
			maxsdev = sdev
			ThresholdWave[1] = {k1}
			//print sdev, k1
		endif
	endfor
end

function FindThresh2(Tab)
wave Tab
	wave ThresholdWave
	variable k1, k2, sdev, maxsdev = 0
	for( k1 = 1; k1 < kNumBins-2; k1 += 1 )
		for( k2 = k1+1; k2 < kNumBins-1; k2 += 1 )
			sdev = Tab[1][k1] + Tab[k1+1][k2] + Tab[k2+1][kNumBins-1]
			if( sdev > maxsdev )
				maxsdev = sdev
				ThresholdWave[1] = {k1, k2}
				//print sdev, k1, k2
			endif
		endfor
	endfor
end

function FindThresh3(Tab)
wave Tab
	wave ThresholdWave
	variable k1, k2, k3, sdev, maxsdev = 0
	for( k1 = 1; k1 < kNumBins-3; k1 += 1 )
		for( k2 = k1+1; k2 < kNumBins-2; k2 += 1 )
			for( k3 = k2+1; k3 < kNumBins-1; k3 += 1 )
				sdev = Tab[1][k1] + Tab[k1+1][k2] + Tab[k2+1][k3] + Tab[k3+1][kNumBins-1]
				if( sdev > maxsdev )
					maxsdev = sdev
					ThresholdWave[1] = {k1, k2, k3}
					//print sdev, k1, k2, k3
				endif
			endfor
		endfor
	endfor
end

function FindThresh4(Tab)
wave Tab
	wave ThresholdWave
	variable k1, k2, k3, k4, sdev, maxsdev = 0
	for( k1 = 1; k1 < kNumBins-4; k1 += 1 )
		for( k2 = k1+1; k2 < kNumBins-3; k2 += 1 )
			for( k3 = k2+1; k3 < kNumBins-2; k3 += 1 )
				for( k4 = k3+1; k4 < kNumBins-1; k4 += 1 )
					sdev = Tab[1][k1] + Tab[k1+1][k2] + Tab[k2+1][k3] + Tab[k3+1][k4] + Tab[k4+1][kNumBins-1]
					if( sdev > maxsdev )
						maxsdev = sdev
						ThresholdWave[1] = {k1, k2, k3, k4}
						//print sdev, k1, k2, k3, k4
					endif
				endfor
			endfor
		endfor
	endfor
end