CCSegStatResultDisplay

#!/usr/bin/env python3

import os
import sys
import getopt
import errno
import csv

import h5py

import numpy
import scipy.stats

import matplotlib.pyplot as plt

import cv2


class CCSegStatResultsDisplayer:
	
	def crossSimple(self, a, b):
		c = numpy.zeros((3, 1))

		c[0] = b.flat[2] * a.flat[1] - b.flat[1] * a.flat[2]
		c[1] = b.flat[0] * a.flat[2] - b.flat[2] * a.flat[0]
		c[2] = b.flat[1] * a.flat[0] - b.flat[0] * a.flat[1]
		
		return c
	
	def tubeSurface(self, vertices, valuesAtVertices, tubeRadius = 50):
		
		a = numpy.take(vertices, [1], axis = 1) - numpy.take(vertices, [0], axis = 1)
		b = numpy.array([[0], [0], [1]])

		c = self.crossSimple(a, b)
		
		if numpy.count_nonzero(c) > 0:
			b = numpy.array([[1], [0], [0]])
			c = self.crossSimple(a, b)
		b = self.crossSimple(c, a)

		normb = numpy.sqrt(numpy.sum(b * b))
		if normb > 0:
			b = b / normb
		
		unitNormals = numpy.zeros_like(vertices)
		unitNormals[:, 0] = numpy.array(b.ravel())
		
		for z in range(vertices.shape[1] - 1):
			a = numpy.take(vertices, [z + 1], axis = 1) - numpy.take(vertices, [z], axis = 1)
			c = self.crossSimple(a, b)
			oldB = numpy.array(b)
			b = self.crossSimple(c, a)

			normb = numpy.sqrt(numpy.sum(b * b))
			if normb > 0:
				b = b / normb
			#print b.T
			unitNormals[:, z + 1] = numpy.array(b.ravel())
#
#		% determine the binormals, which is B = T x N
#		% T = tangent
#		% N = normal
		tangentVectors = vertices - numpy.roll(vertices, 1, axis = 1)
		
		binormalVectors = numpy.zeros_like(vertices)
		binormalVectors[0] = tangentVectors[1] * unitNormals[2] - tangentVectors[2] * unitNormals[1] 
		binormalVectors[1] = tangentVectors[2] * unitNormals[0] - tangentVectors[0] * unitNormals[2] 
		binormalVectors[2] = tangentVectors[0] * unitNormals[1] - tangentVectors[1] * unitNormals[0] 

		speed = numpy.sqrt(numpy.sum(tangentVectors * tangentVectors, axis = 0))
		
		binormalVectors = binormalVectors / speed
		
		print(binormalVectors)


		quit()		
#
#		BiNormalVectors = bsxfun(@rdivide, cross(TangentVectors, unitnormals, 2), speed);
#		% hold off;
#		% plot3(XYZ(:, 1), XYZ(:, 2), XYZ(:, 3), '*-');
#		% hold on;
#		% h = quiver3(XYZ(:, 1), XYZ(:, 2), XYZ(:, 3), unitnormals(:, 1), unitnormals(:, 2), unitnormals(:, 3));
#		% set(h, 'Color', 'r');
#		% h = quiver3(XYZ(:, 1), XYZ(:, 2), XYZ(:, 3), BiNormalVectors(:, 1), BiNormalVectors(:, 2), BiNormalVectors(:, 3));
#		% set(h, 'Color', 'b');
#
#		N = 20;
#		ANGLES = linspace(0, 2 * pi, N + 1);
#		ANGLES = ANGLES(1:N)';
#
#		VCell = cell(size(XYZ, 1), 1);
#		FCell = cell(size(XYZ, 1), 1);
#
#		COSANGLES = cos(ANGLES);
#		SINANGLES = sin(ANGLES);
#		for z = 1:size(XYZ, 1)
#			VR = (bsxfun(@times, COSANGLES, unitnormals(z, :)) + bsxfun(@times, SINANGLES, BiNormalVectors(z, :))) * R / 2;
#			VCell{z} = bsxfun(@plus, XYZ(z, :), VR);
#			
#			I = (1:N)';
#			% for each vertex make faces
#			MinusOneLevel = I - 1;
#			MinusOneLevel(MinusOneLevel == 0) = N;
#			%F = 
#			FCell{z} = ...
#				[I + (z - 1) * N, ... same level, same column
#				 mod(I, N) + 1 + (z - 1) * N, ... one level up, same column
#				 mod(I + z * N - 1, size(XYZ, 1) * N) + 1; ... same level, next column
#				 I + (z - 1) * N, ... same level, same column
#				 mod(I + z * N - 1, size(XYZ, 1) * N) + 1, ... same level, next column
#				 mod(MinusOneLevel + z * N - 1, size(XYZ, 1) * N) + 1]; ... one level down, next column
#				
#			%FCell{z} = ...
#			%	[I, mod((I + N) - 1, N * size(XYZ, 1)) + 1, mod(mod(I, N) + z * N, N * size(XYZ, 1)) + 1];
#			%keyboard;
#			% first row, I, I - 1, I + N
#			% second row, I, I + N, I + N - 1
#		end
#
#		%streamline(VCell);
#		F = cat(1, FCell{:});
#		V = cat(1, VCell{:});
#
#		%II = [repmat(1, N, 1); repmat(2, N, 1)];
#		%CData = repmat(II(:)', size(XYZ, 1), 1);
#		%CData = repmat(II(:), 1, size(XYZ, 1));
#		CData = repmat(Data(:)', N, 1);
#		CData = CData(:);
#
#		% patch('Faces', F, 'Vertices', V, 'FaceColor', 'r', 'EdgeColor', 'k');
#		% axis equal;
#		% keyboard;
#
#
#		function c=crossSimple(a,b)
#
#		c(1) = b(3)*a(2) - b(2)*a(3);
#		c(2) = b(1)*a(3) - b(3)*a(1);
#		c(3) = b(2)*a(1) - b(1)*a(2);	

	def plotData(self, data, plotType = 'pvalue', groupLabels = None, legendText = None):

		numNodes = numpy.size(data)
		
		interpRange = numpy.array([numpy.min(self.idealCCInterpPoints[self.maskClosed]), numpy.max(self.idealCCInterpPoints[self.maskClosed])])
		
		self.idealCCInterpPoints[self.maskClosed] = 1.0 - (self.idealCCInterpPoints[self.maskClosed] - interpRange[0]) / (interpRange[1] - interpRange[0])
		self.idealCCInterpPoints[self.maskClosed] = self.idealCCInterpPoints[self.maskClosed] * (numNodes - 1.0) + 1.0

		self.tubeSurface(numpy.concatenate((self.boundaryContour, numpy.zeros((1, self.boundaryContour.shape[1]))), axis = 0), data)
	
	def __init__(self):

		#TODO: put in non-hardcoded location for this file, DONE
		#idealCCFile = 'ideal_cc.hdf5'
		scriptPath = os.path.realpath(__file__)
		(head, tail) = os.path.split(scriptPath)

		idealCCFile = os.path.join(head, 'data', 'ideal_cc.hdf5')

		if not os.path.isfile(idealCCFile):
			print("ideal CC file not found")
			quit()
		
		FID = h5py.File(idealCCFile, 'r')
		self.idealCCInterpPoints = numpy.array(FID["InterpPoints"]).T
		FID.close()	
		
		self.maskClosed = numpy.nonzero(self.idealCCInterpPoints)
		

		#print interpRange

		T = numpy.uint8(self.idealCCInterpPoints > 0)
		
		T = numpy.pad(T, [1, 1], mode='constant', constant_values=0)
		#print T 
		#print T.dtype
		boundaryContour, hierarchy = CCSegUtils.findContours(T, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
		del hierarchy
		del T

		# the minus one is due to the padding we added
		self.boundaryContour = numpy.squeeze(boundaryContour[0].T) - 1
		self.boundaryContour = self.boundaryContour[:, 0::5]
		#print boundaryContour
		#print boundaryContour.shape
		
		#plt.imshow(T)
		
		#plt.colorbar()
		#plt.gcf().set_size_inches((20, 10), forward = True)
		#plt.show()
		
		#plt.show()
		# set up

		#self.plotType = plotType

		#quit()




def usage():
	print("Usage: " + sys.argv[0] + " [options] <stattest hdf5 file>")
	print()
	print("\tDESCRIPTION")
	print()
	print("\t\tDisplays statistical analysis results on callosal thickness profiles created by CCSegStatTest")
	print()
	print("\tARGUMENTS")
	print()
	print("\t\t<stattest hdf5 file>: the hdf5 file produced by CCSegStatTest, must contain the following variables")
	print("\t\t\tGroupLabels")
	print("\t\t\tGroupIDX")
	print("\t\t\tobservedP")
	print("\t\t\tobservedT")
	print("\t\t\tpermP")
	print("\t\t\tomnibusP")
	print("\t\t\tFDRP")

def main():
	opts, args = getopt.getopt(sys.argv[1:], "h", [])

	numpy.set_printoptions(precision = 3, formatter = {'float':lambda x: "%.3f" % x})

	inputHDF5File = None

	if len(args) != 1:
		print("The number of arguments must be 1")
		usage()
		exit()

	inputHDF5File = args[0]

	if not os.path.isfile(inputHDF5File):
		print("The input file does not exist")
		quit()

	#for o, a in opts:
	#	if o == '--cull-percent':
	#		try:
	#			cullPercent = float(a)
	#		except Exception:
	#			print "Cull percentage was not formatted in a valid way, it must be a floating point number"
	#			quit()

	statOutput = dict()

	FID = h5py.File(inputHDF5File, 'r')

	statOutput["permP"] = numpy.array(FID["permP"])
	statOutput["FDRP"] = numpy.array(FID["FDRP"])
	statOutput["observedP"] = numpy.array(FID["observedP"])
	statOutput["observedT"] = numpy.array(FID["observedT"])
	statOutput["omnibusP"] = numpy.array(FID["omnibusP"])
	statOutput["groupLabels"] = numpy.array(FID['groupLabels']).tolist()
	FID.close()
	
	numNodes = numpy.size(statOutput["permP"])

	resultsDisplayer = CCSegStatResultsDisplayer()
	resultsDisplayer.plotData(data = statOutput["permP"])

	print(sys.argv[0])

if __name__ == "__main__":
	main()