First commit

AddLayers_To3DVol.py

+#!/usr/bin/env python3
+
+
+import imageio 
+import numpy as np
+
+import sys
+import argparse
+import time 
+
+def arg_parser():
+	parser = argparse.ArgumentParser(description='Concatenating layer to 3D volume')
+	required = parser.add_argument_group('Required')
+	required.add_argument('--input', type=str, required=True,
+						  help='3D TIFF file (multi-layer)')
+	required.add_argument('--layer1', type=str, required=True,
+						  help='2D TIFF file (single layer)')
+	required.add_argument('--output', type=str, required=True,
+						  help='Combined TIFF file (multi-layer)')
+	options = parser.add_argument_group('Options')
+	required.add_argument('--layer2', type=str, required=True,
+						  help='2D TIFF file (single layer)')
+	options.add_argument('--repeat', type=int, default=15,
+						  help='Repeat tile number (default 15)')
+	return parser
+
+
+#MAIN
+def main(args=None):
+	args = arg_parser().parse_args(args)
+
+	Time1 = time.time()
+	#print(bcolors.BOLDWHITE+"Time1: "+str(Time1)+bcolors.ENDC)
+
+	imgInput_name = args.input
+	imgLayer1_name = args.layer1
+	imgLayer2_name = args.layer2
+	output_name = args.output
+	repeatNb = args.repeat
+	
+	# Read Combined image - all layers
+	print("\nReading input file - 3D volume...")
+	imgInput = imageio.mimread(imgInput_name,memtest=False)
+	imgInput = np.array(imgInput)
+	# print('\nimgInput type: ', imgInput.dtype)
+	print('\t imgInput shape: ', imgInput.shape)
+
+	# Read Layer1 image
+	print("Reading Layer1 file...")
+	imgLayer1 = imageio.imread(imgLayer1_name)
+	# print('imgLayer1 type: ', imgLayer1.dtype)
+	print('\t imgLayer1 shape: ', imgLayer1.shape)
+
+	# Resample layer1 (repeating values, to match input size)
+	imgLayer1_repeat0 = np.repeat(imgLayer1, repeatNb, axis=0)
+	imgLayer1_repeat = np.repeat(imgLayer1_repeat0, repeatNb, axis=1)
+	imgLayer1_repeat = np.expand_dims(imgLayer1_repeat, axis=0)
+	print('\t imgLayer1_repeat shape: ', imgLayer1_repeat.shape)
+
+	# Stack layer to imgInput, to generate one 3D volume
+	print("Adding 2D layer to 3D volume...")
+	imgAll = np.concatenate((imgInput,imgLayer1_repeat), axis = 0)
+
+	if args.layer2:
+		print("Reading Layer2 file...")
+		imgLayer2 = imageio.imread(imgLayer2_name)
+		# print('imgLayer2 type: ', imgLayer2.dtype)
+		print('\t imgLayer2 shape: ', imgLayer2.shape)
+
+		imgLayer2_repeat0 = np.repeat(imgLayer2, repeatNb, axis=0)
+		imgLayer2_repeat = np.repeat(imgLayer2_repeat0, repeatNb, axis=1)
+		imgLayer2_repeat = np.expand_dims(imgLayer2_repeat, axis=0)
+		print('\t imgLayer2_repeat shape: ', imgLayer2_repeat.shape)
+
+		# Stack layer to imgInput, to generate one 3D volume
+		print("Adding 2D layer to 3D volume...")
+		imgAll = np.concatenate((imgAll,imgLayer2_repeat), axis = 0)
+
+	print("Saving output volume...")
+	print('\t imgAll shape: ', imgAll.shape)
+	imageio.mimwrite(output_name,imgAll)
+
+	Time2 = time.time()
+	TimeDiff = Time2 - Time1
+	print("Execution Time: "+str(TimeDiff))
+
+if __name__ == '__main__':
+	sys.exit(main(sys.argv[1:]))
+

Compute_Density.py

+#!/usr/bin/env python3
+
+
+import imageio 
+import numpy as np
+import pandas as pd
+from sklearn.metrics import mean_squared_error
+
+import sys
+import argparse
+import time 
+import os 
+
+def arg_parser():
+	parser = argparse.ArgumentParser(description='Data analysis - density computation')
+	required = parser.add_argument_group('Required')
+	required.add_argument('--pred', type=str, required=True,
+		help='Prediction TIFF file (single layer)')
+	required.add_argument('--groundtruth', type=str, required=True,
+		help='Ground Truth TIFF file (single layer)')
+	required.add_argument('--adjtilesdim', type=int, required=True,
+		help='Adjacent tiles dimensions (e.g. 1, 3 or 5) to exclude NaN border')
+	required.add_argument('--output', type=str, required=True,
+		help='CSV file')
+	options = parser.add_argument_group('Options')
+	options.add_argument('--inclusionmask', type=str, 
+						  help='save inclusion mask as output (diff < 2 pixels)')
+	options.add_argument('--exclusionmask', type=str, 
+						  help='save exclusion mask as output (diff >= 2 pixels)')
+	options.add_argument('--threshold', type=float, default = 2.0,
+		help='threshold for inclusion / exclusion mask')
+	options.add_argument('--verbose', action="store_true",
+						  help='verbose mode')
+	return parser
+
+#MAIN
+def main(args=None):
+	args = arg_parser().parse_args(args)
+	# print(args)
+
+	Time1 = time.time()
+	#print(bcolors.BOLDWHITE+"Time1: "+str(Time1)+bcolors.ENDC)
+
+	imgPred_name = args.pred
+	imgPred_basename = os.path.basename(imgPred_name)
+	AdjacentTilesDim = args.adjtilesdim
+	imgGT_name = args.groundtruth
+	output_name = args.output
+	inclusion_mask_name = args.inclusionmask
+	exclusion_mask_name = args.exclusionmask
+	threshold = args.threshold
+	
+
+	# Read Pred and GroundTruth images
+	#print("Reading Pred file...")
+	imgPred = imageio.imread(imgPred_name)
+
+
+	#print("Reading GroundTruth file...")
+	imgGT = imageio.imread(imgGT_name)
+
+	# Remove NaN border when needed	
+	imgPred_Crop = imgPred
+	imgGT_Crop = imgGT
+	if (AdjacentTilesDim == 3):
+		Border = 1
+		imgPred_Crop = imgPred[Border:-Border,Border:-Border]
+		imgGT_Crop = imgGT[Border:-Border,Border:-Border]
+	elif (AdjacentTilesDim == 5):
+		Border = 2
+		imgPred_Crop = imgPred[Border:-Border,Border:-Border]
+		imgGT_Crop = imgGT[Border:-Border,Border:-Border]
+
+	# Quality control - NaN
+	TestNaN_imgPred_Crop = np.any(np.isnan(imgPred_Crop))
+	TestNaN_imgGT_Crop = np.any(np.isnan(imgGT_Crop))
+
+	# #  List indices with Nan values
+	# ListNaN_imgPred_Crop = np.argwhere(np.isnan(imgPred_Crop))
+	# print('ListNaN_imgPred_Crop: ', ListNaN_imgPred_Crop)
+
+	# Verbose mode
+	if args.verbose:
+		print('imgPred type: ', imgPred.dtype)
+		print('imgPred shape: ', imgPred.shape)		
+		print('imgGT type: ', imgGT.dtype)
+		print('imgGT shape: ', imgGT.shape)
+		print('imgPred_Crop shape: ', imgPred_Crop.shape)
+		print('imgGT_Crop shape: ', imgGT_Crop.shape)
+		print('TestNaN_imgPred_Crop: ', TestNaN_imgPred_Crop)
+		print('TestNaN_imgGT_Crop: ', TestNaN_imgGT_Crop)
+
+	# Compute Image difference
+	np_diff = np.abs(imgPred_Crop - imgGT_Crop)
+	print('np_diff type: ', np_diff.dtype)
+	print('np_diff shape: ', np_diff.shape)		
+		
+	# Generate output mask
+	np_exclusionmask = np.uint8(np.where(np_diff >= threshold, 1, 0))
+	np_inclusionmask = np.uint8(np.where(np_diff < threshold, 1, 0))
+	print('np_inclusionmask type: ', np_inclusionmask.dtype)
+	print('np_inclusionmask shape: ', np_inclusionmask.shape)		
+	
+	# Compute density
+	density = np.sum(np_inclusionmask) / (np_inclusionmask.shape[0] * np_inclusionmask.shape[1])
+	print('density: ',density)
+
+	QC_data = np.array([[imgPred_basename,density]])
+	Columns = ['FileName','Density']
+	df = pd.DataFrame(QC_data,columns=Columns)
+	print(df.head())
+	df.to_csv(output_name, index=False)
+
+	if inclusion_mask_name is not None:
+		print('\t\t Writing output inclusion mask - imageio...')
+		imageio.imwrite(inclusion_mask_name, np_inclusionmask)
+
+	if inclusion_mask_name is not None:
+		print('\t\t Writing output exclusion mask - imageio...')
+		imageio.imwrite(exclusion_mask_name, np_exclusionmask)
+
+	Time2 = time.time()
+	#print(bcolors.BOLDWHITE+"Time2: "+str(Time2)+bcolors.ENDC)
+	TimeDiff = Time2 - Time1
+	#print("Computing Time: "+str(TimeDiff))
+
+if __name__ == '__main__':
+	sys.exit(main(sys.argv[1:]))
+

Compute_DispLMA_Mask.py

+#!/usr/bin/env python3
+
+
+import imageio 
+import numpy as np
+import pandas as pd
+
+import sys
+import argparse
+import time 
+import os 
+
+def arg_parser():
+	parser = argparse.ArgumentParser(description='Data analysis - NaN computation')
+	required = parser.add_argument_group('Required')
+	required.add_argument('--disp_lma', type=str, required=True,
+		help='Prediction TIFF file (single layer)')
+	required.add_argument('--output', type=str, required=True,
+		help='Output CSV file')
+	options = parser.add_argument_group('Options')
+	options.add_argument('--verbose', action="store_true",
+						  help='verbose mode')
+	options.add_argument('--mask', type=str,
+						  help='verbose mode')
+	return parser
+
+
+#MAIN
+def main(args=None):
+	args = arg_parser().parse_args(args)
+
+	Time1 = time.time()
+	#print(bcolors.BOLDWHITE+"Time1: "+str(Time1)+bcolors.ENDC)
+
+	DispLMA_name = args.disp_lma
+	output_name = args.output
+
+	# Read Pred and GroundTruth images
+	#print("Reading Pred file...")
+	imgDispLMA = imageio.imread(DispLMA_name)
+
+	# Verbose mode
+	if args.verbose:
+		print('imgDispLMA type: ', imgDispLMA.dtype)
+		print('imgDispLMA shape: ', imgDispLMA.shape)		
+		# print('imgGT type: ', imgGT.dtype)
+
+	# Compute NaN
+	Bool_NaN = np.isnan(imgDispLMA)
+	Nb_NaN = np.count_nonzero(Bool_NaN)
+	Nb_NotNaN = np.count_nonzero(~Bool_NaN)
+	Nb_Total = imgDispLMA.shape[0] * imgDispLMA.shape[1]
+	NaN_Percent = Nb_NaN / Nb_Total
+	
+	if args.verbose:
+		print('Nb_NaN: ', Nb_NaN)
+		print('Nb_NotNaN: ', Nb_NotNaN)
+		print('Nb_Total: ', Nb_Total)
+		print('NaN_Percent: ', NaN_Percent)
+
+
+	# Mask
+	DispLMA_Mask = np.uint8(np.reshape(~Bool_NaN, imgDispLMA.shape))
+	if args.verbose:
+		print('DispLMA_Mask type: ', DispLMA_Mask.dtype)
+		print('DispLMA_Mask shape: ', DispLMA_Mask.shape)		
+	
+
+	QC_data = np.array([[DispLMA_name, Nb_NaN, Nb_NotNaN, Nb_Total, NaN_Percent]])
+	Columns = ['FileName','Nb_NaN', 'Nb_NotNaN', 'Nb_Total', 'NaN_Percent']
+	df = pd.DataFrame(QC_data,columns=Columns)
+	print(df.head())
+	df.to_csv(output_name, index=False)
+
+	if (args.mask is not None):
+		Mask_FileName = args.mask
+		imageio.imwrite(Mask_FileName, DispLMA_Mask * 255)
+
+
+	Time2 = time.time()
+	#print(bcolors.BOLDWHITE+"Time2: "+str(Time2)+bcolors.ENDC)
+	TimeDiff = Time2 - Time1
+	#print("Computing Time: "+str(TimeDiff))
+
+if __name__ == '__main__':
+	sys.exit(main(sys.argv[1:]))
+

Compute_RMSE.py

+#!/usr/bin/env python3
+
+
+import imageio 
+import numpy as np
+import pandas as pd
+from sklearn.metrics import mean_squared_error
+
+import sys
+import argparse
+import time 
+import os 
+
+def arg_parser():
+	parser = argparse.ArgumentParser(description='Data analysis - RMSE computation')
+	required = parser.add_argument_group('Required')
+	required.add_argument('--pred', type=str, required=True,
+		help='Prediction TIFF file (single layer)')
+	required.add_argument('--groundtruth', type=str, required=True,
+		help='Ground Truth TIFF file (single layer)')
+	required.add_argument('--adjtilesdim', type=int, required=True,
+		help='Adjacent tiles dimensions (e.g. 1, 3 or 5) to exclude NaN border')
+	required.add_argument('--output', type=str, required=True,
+		help='CSV file')
+	options = parser.add_argument_group('Options')
+	options.add_argument('--verbose', action="store_true",
+						  help='verbose mode')
+	return parser
+
+#MAIN
+def main(args=None):
+	args = arg_parser().parse_args(args)
+
+	Time1 = time.time()
+	#print(bcolors.BOLDWHITE+"Time1: "+str(Time1)+bcolors.ENDC)
+
+	imgPred_name = args.pred
+	imgPred_basename = os.path.basename(imgPred_name)
+	AdjacentTilesDim = args.adjtilesdim
+	imgGT_name = args.groundtruth
+	output_name = args.output
+
+	# Read Pred and GroundTruth images
+	#print("Reading Pred file...")
+	imgPred = imageio.imread(imgPred_name)
+
+
+	#print("Reading GroundTruth file...")
+	imgGT = imageio.imread(imgGT_name)
+
+	# Remove NaN border when needed	
+	imgPred_Crop = imgPred
+	imgGT_Crop = imgGT
+	if (AdjacentTilesDim == 3):
+		Border = 1
+		imgPred_Crop = imgPred[Border:-Border,Border:-Border]
+		imgGT_Crop = imgGT[Border:-Border,Border:-Border]
+	elif (AdjacentTilesDim == 5):
+		Border = 2
+		imgPred_Crop = imgPred[Border:-Border,Border:-Border]
+		imgGT_Crop = imgGT[Border:-Border,Border:-Border]
+
+	# Quality control - NaN
+	TestNaN_imgPred_Crop = np.any(np.isnan(imgPred_Crop))
+	TestNaN_imgGT_Crop = np.any(np.isnan(imgGT_Crop))
+
+	# #  List indices with Nan values
+	# ListNaN_imgPred_Crop = np.argwhere(np.isnan(imgPred_Crop))
+	# print('ListNaN_imgPred_Crop: ', ListNaN_imgPred_Crop)
+
+	# Verbose mode
+	if args.verbose:
+		# print('imgPred type: ', imgPred.dtype)
+		print('imgPred shape: ', imgPred.shape)		
+		# print('imgGT type: ', imgGT.dtype)
+		print('imgGT shape: ', imgGT.shape)
+		print('imgPred_Crop shape: ', imgPred_Crop.shape)
+		print('imgGT_Crop shape: ', imgGT_Crop.shape)
+		print('TestNaN_imgPred_Crop: ', TestNaN_imgPred_Crop)
+		print('TestNaN_imgGT_Crop: ', TestNaN_imgGT_Crop)
+
+	# Compute RMSE
+	rmse = mean_squared_error(imgGT_Crop, imgPred_Crop, squared=False)
+	print('rmse: ',rmse)
+
+	QC_data = np.array([[imgPred_basename,rmse]])
+	Columns = ['FileName','RMSE']
+	df = pd.DataFrame(QC_data,columns=Columns)
+	print(df.head())
+	df.to_csv(output_name, index=False)
+
+	Time2 = time.time()
+	#print(bcolors.BOLDWHITE+"Time2: "+str(Time2)+bcolors.ENDC)
+	TimeDiff = Time2 - Time1
+	#print("Computing Time: "+str(TimeDiff))
+
+if __name__ == '__main__':
+	sys.exit(main(sys.argv[1:]))
+

Compute_RMSE_WithFiltering.py

+#!/usr/bin/env python3
+
+
+import imageio 
+import numpy as np
+import pandas as pd
+from sklearn.metrics import mean_squared_error
+
+import sys
+import argparse
+import time 
+import os 
+
+def arg_parser():
+	parser = argparse.ArgumentParser(description='Data analysis - RMSE computation')
+	required = parser.add_argument_group('Required')
+	required.add_argument('--pred', type=str, required=True,
+		help='Prediction TIFF file (single layer)')
+	required.add_argument('--groundtruth', type=str, required=True,
+		help='Ground Truth TIFF file (single layer)')
+	required.add_argument('--confidence', type=str, required=True,
+		help='Confidence TIFF file (single layer)')
+	required.add_argument('--disp_lma', type=str, required=True,
+		help='LMA disparity TIFF file (single layer)')
+	required.add_argument('--adjtilesdim', type=int, required=True,
+		help='Adjacent tiles dimensions (e.g. 1, 3 or 5) to exclude NaN border')
+	required.add_argument('--output', type=str, required=True,
+		help='CSV file')
+	options = parser.add_argument_group('Options')
+	options.add_argument('--threshold', type=float, default = 0.15,
+						  help='threshold on confidence map')
+	options.add_argument('--output_mask', type=str,
+		help='output mask image (TIFF file)')
+	options.add_argument('--verbose', action="store_true",
+						  help='verbose mode')
+	return parser
+
+# Remove NaN border
+def cropping(img, border):
+	if border == 0:
+		return img
+	else:
+		img_cropped = img[border:-border,border:-border]
+		return img_cropped
+
+#MAIN
+def main(args=None):
+	args = arg_parser().parse_args(args)
+	
+	Time1 = time.time()
+	#print(bcolors.BOLDWHITE+"Time1: "+str(Time1)+bcolors.ENDC)
+
+	imgPred_name = args.pred
+	imgPred_basename = os.path.basename(imgPred_name)
+	AdjacentTilesDim = args.adjtilesdim
+	imgGT_name = args.groundtruth
+	imgDispLMA_name = args.disp_lma
+	imgConfidence_name = args.confidence
+	output_name = args.output
+	outputmask_name = args.output_mask
+	threshold = args.threshold
+
+	# Read images
+	#print("Reading Pred file...")
+	imgPred = imageio.imread(imgPred_name)
+
+	#print("Reading GroundTruth file...")
+	imgGT = imageio.imread(imgGT_name)
+
+	#print("Reading Confidence file...")
+	imgConfidence = imageio.imread(imgConfidence_name)
+
+	#print("Reading DispLMA file...")
+	imgDispLMA = imageio.imread(imgDispLMA_name)
+
+	# Remove NaN border when needed	
+	if (AdjacentTilesDim == 3):
+		Border = 1	
+	elif (AdjacentTilesDim == 5):
+		Border = 2
+	else:
+		Border = 0
+	imgPred_Crop = cropping(imgPred, Border)
+	imgGT_Crop = cropping(imgGT, Border)
+	imgConfidence_Crop = cropping(imgConfidence, Border)
+	imgDispLMA_Crop = cropping(imgDispLMA, Border)
+
+	# Quality control - NaN
+	TestNaN_imgPred_Crop = np.any(np.isnan(imgPred_Crop))
+	TestNaN_imgGT_Crop = np.any(np.isnan(imgGT_Crop))
+
+	# #  List indices with Nan values
+	# ListNaN_imgPred_Crop = np.argwhere(np.isnan(imgPred_Crop))
+	# print('ListNaN_imgPred_Crop: ', ListNaN_imgPred_Crop)
+
+	# Verbose mode
+	if args.verbose:
+		# print('imgPred type: ', imgPred.dtype)
+		print('imgPred shape: ', imgPred.shape)		
+		# print('imgGT type: ', imgGT.dtype)
+		print('imgGT shape: ', imgGT.shape)
+		print('imgPred_Crop shape: ', imgPred_Crop.shape)
+		print('imgGT_Crop shape: ', imgGT_Crop.shape)
+		print('TestNaN_imgPred_Crop: ', TestNaN_imgPred_Crop)
+		print('TestNaN_imgGT_Crop: ', TestNaN_imgGT_Crop)
+
+
+	# Define sample_weight using Confidence and DispLMA maps
+	imgDispLMAMask_Crop = np.uint8(np.reshape(~np.isnan(imgDispLMA_Crop), imgDispLMA_Crop.shape))
+	imgConfidenceMask_Crop = np.uint8(np.where(imgConfidence_Crop >= threshold, 1, 0))
+	imgSampleWeight_Crop = np.uint8(np.logical_and(imgDispLMAMask_Crop, imgConfidenceMask_Crop))
+
+	if outputmask_name is not None:
+		imageio.imwrite(outputmask_name, imgSampleWeight_Crop * 255)
+
+	# Compute RMSE
+	rmse = mean_squared_error(imgGT_Crop, imgPred_Crop, sample_weight=imgSampleWeight_Crop, squared=False)
+	print('rmse: ',rmse)
+
+	QC_data = np.array([[imgPred_basename,rmse]])
+	Columns = ['FileName','RMSE']
+	df = pd.DataFrame(QC_data,columns=Columns)
+	#print(df.head())
+
+	print('Saving CSV file...')
+	df.to_csv(output_name, index=False)
+
+	Time2 = time.time()
+	#print(bcolors.BOLDWHITE+"Time2: "+str(Time2)+bcolors.ENDC)
+	TimeDiff = Time2 - Time1
+	#print("Computing Time: "+str(TimeDiff))
+
+if __name__ == '__main__':
+	sys.exit(main(sys.argv[1:]))
+

Compute_RMSE_WithThreshold.py

+#!/usr/bin/env python3
+
+
+import imageio 
+import numpy as np
+import pandas as pd
+from sklearn.metrics import mean_squared_error
+
+import sys
+import argparse
+import time 
+import os 
+
+def arg_parser():
+	parser = argparse.ArgumentParser(description='Data analysis - RMSE computation')
+	required = parser.add_argument_group('Required')
+	required.add_argument('--pred', type=str, required=True,
+		help='Prediction TIFF file (single layer)')
+	required.add_argument('--groundtruth', type=str, required=True,
+		help='Ground Truth TIFF file (single layer)')
+	required.add_argument('--adjtilesdim', type=int, required=True,
+		help='Adjacent tiles dimensions (e.g. 1, 3 or 5) to exclude NaN border')
+	required.add_argument('--output', type=str, required=True,
+		help='CSV file')
+	options = parser.add_argument_group('Options')
+	options.add_argument('--threshold', type=float, default = 2.0,
+						  help='threshold on image difference')
+	options.add_argument('--output_mask', type=str,
+		help='output mask image (TIFF file)')
+	options.add_argument('--verbose', action="store_true",
+						  help='verbose mode')
+	return parser
+
+# Remove NaN border
+def cropping(img, border):
+	if border == 0:
+		return img
+	else:
+		img_cropped = img[border:-border,border:-border]
+		return img_cropped
+
+#MAIN
+def main(args=None):
+	args = arg_parser().parse_args(args)
+	
+	Time1 = time.time()
+	#print(bcolors.BOLDWHITE+"Time1: "+str(Time1)+bcolors.ENDC)
+
+	imgPred_name = args.pred
+	imgPred_basename = os.path.basename(imgPred_name)
+	AdjacentTilesDim = args.adjtilesdim
+	imgGT_name = args.groundtruth
+	output_name = args.output
+	outputmask_name = args.output_mask
+	threshold = args.threshold
+
+	# Read images
+	#print("Reading Pred file...")
+	imgPred = imageio.imread(imgPred_name)
+
+	#print("Reading GroundTruth file...")
+	imgGT = imageio.imread(imgGT_name)
+
+	# Remove NaN border when needed	
+	if (AdjacentTilesDim == 3):
+		Border = 1	
+	elif (AdjacentTilesDim == 5):
+		Border = 2
+	else:
+		Border = 0
+	imgPred_Crop = cropping(imgPred, Border)
+	imgGT_Crop = cropping(imgGT, Border)
+	
+	# Quality control - NaN
+	TestNaN_imgPred_Crop = np.any(np.isnan(imgPred_Crop))
+	TestNaN_imgGT_Crop = np.any(np.isnan(imgGT_Crop))
+
+	# #  List indices with Nan values
+	# ListNaN_imgPred_Crop = np.argwhere(np.isnan(imgPred_Crop))
+	# print('ListNaN_imgPred_Crop: ', ListNaN_imgPred_Crop)
+
+	# Verbose mode
+	if args.verbose:
+		# print('imgPred type: ', imgPred.dtype)
+		print('imgPred shape: ', imgPred.shape)		
+		# print('imgGT type: ', imgGT.dtype)
+		print('imgGT shape: ', imgGT.shape)
+		print('imgPred_Crop shape: ', imgPred_Crop.shape)
+		print('imgGT_Crop shape: ', imgGT_Crop.shape)
+		print('TestNaN_imgPred_Crop: ', TestNaN_imgPred_Crop)
+		print('TestNaN_imgGT_Crop: ', TestNaN_imgGT_Crop)
+
+
+	# Define sample_weight, where absolute difference is below threshold (2 pixels), to remove outliers
+	imgAbsDiff_Crop = np.absolute(imgPred_Crop - imgGT_Crop)
+	imgSampleWeight_Crop = np.uint8(np.where(imgAbsDiff_Crop < threshold, 1, 0))
+	
+	if outputmask_name is not None:
+		imageio.imwrite(outputmask_name, imgSampleWeight_Crop * 255)
+
+	# Compute RMSE
+	rmse = mean_squared_error(imgGT_Crop, imgPred_Crop, sample_weight=imgSampleWeight_Crop, squared=False)
+	print('rmse: ',rmse)
+
+	QC_data = np.array([[imgPred_basename,rmse]])
+	Columns = ['FileName','RMSE_WithThreshold']
+	df = pd.DataFrame(QC_data,columns=Columns)
+	#print(df.head())
+
+	print('Saving CSV file...')
+	df.to_csv(output_name, index=False)
+
+	Time2 = time.time()
+	#print(bcolors.BOLDWHITE+"Time2: "+str(Time2)+bcolors.ENDC)
+	TimeDiff = Time2 - Time1
+	#print("Computing Time: "+str(TimeDiff))
+
+if __name__ == '__main__':
+	sys.exit(main(sys.argv[1:]))
+

Preprocessing_2DPhaseCorrelation.py

+#!/usr/bin/env python3
+
+
+
+'''
+  Notes:
+	- Pillow 5.1.0. Version 4.1.1 throws error (VelueError):
+	  ~$ (sudo) pip3 install Pillow --upgrade
+	  ~$ python3
+	  >>> import PIL
+	  >>> PIL.PILLOW_VERSION
+	  '5.1.0'
+'''
+
+from PIL import Image
+import xml.etree.ElementTree as ET
+import numpy as np
+import matplotlib.pyplot as plt
+
+import os
+import sys
+import xml.dom.minidom as minidom
+import time
+
+# Added by Clement
+import imageio
+#import tifffile
+import argparse
+
+#http://stackoverflow.com/questions/287871/print-in-terminal-with-colors-using-python
+class bcolors:
+	HEADER = '\033[95m'
+	OKBLUE = '\033[94m'
+	OKGREEN = '\033[92m'
+	WARNING = '\033[38;5;214m'
+	FAIL = '\033[91m'
+	ENDC = '\033[0m'
+	BOLD = '\033[1m'
+	BOLDWHITE = '\033[1;37m'
+	UNDERLINE = '\033[4m'
+
+# Data cropping (border of specific size)
+def crop_border(image, size):
+	List_FirstRows = list(range(0,size))
+	List_LastRows = list(range(image.shape[0]-size,image.shape[0]))
+	List_Rows = List_FirstRows + List_LastRows
+	#print('ListRows: ', List_Rows)
+
+	List_FirstCols = list(range(0,size))
+	List_LastCols = list(range(image.shape[1]-size,image.shape[1]))
+	List_Cols = List_FirstCols + List_LastCols
+	#print('List_Cols: ', List_Cols)
+
+	CroppedImage = np.delete(image,List_Rows,axis=0)
+	CroppedImage = np.delete(CroppedImage,List_Cols,axis=1)
+
+	return CroppedImage
+
+# 2D Corr data cropping
+def delete_tile_splits(image, tile_size):
+	List_RowSplit = list(range(tile_size-1,image.shape[0],tile_size))
+	List_ColSplit = list(range(tile_size-1,image.shape[1],tile_size))
+	#print('List_RowSplit: ', List_RowSplit)
+	#print('List_ColSplit: ', List_ColSplit)
+
+	CroppedImage = np.delete(image,List_RowSplit,axis=0)
+	CroppedImage = np.delete(CroppedImage,List_ColSplit,axis=1)
+
+	return CroppedImage
+
+def save_metadata_layer(image, filename, index, is_cropping, verbose):
+		#metadata_layer = tiles_meta[0,:,:,index]
+		metadata_layer = image.getMeta(index = index, layer = 0)
+		if is_cropping:
+			metadata_layer = crop_border(metadata_layer, 1)
+		if verbose:
+			print('\metadata_layer shape: ', metadata_layer.shape)
+		if not os.path.exists(os.path.dirname(filename)):
+			os.makedirs(os.path.dirname(filename))
+		imageio.imwrite(filename, metadata_layer)
+
+# reshape to tiles (e.g. 16x16)
+def get_tile_images(image, width=8, height=8):
+	_nrows, _ncols, depth = image.shape
+	_size = image.size
+	_strides = image.strides
+	
+	nrows, _m = divmod(_nrows, height)
+	ncols, _n = divmod(_ncols, width)
+	if _m != 0 or _n != 0:
+		return None
+
+	return np.lib.stride_tricks.as_strided(
+		np.ravel(image),
+		shape=(nrows, ncols, height, width, depth),
+		strides=(height * _strides[0], width * _strides[1], *_strides),
+		writeable=False
+	)
+
+# TiffFile has no len exception
+#import imageio
+
+#from libtiff import TIFF
+'''
+Description:
+	Reads TIFF files with multiple layers that were saved by imagej
+Methods:
+	.getstack(items=[])
+		returns np.array, layers are stacked along depth - think of RGB channels
+		@items - if empty = all, if not - items[i] - can be layer index or layer's label name
+	.channel(index)
+		returns np.array of a single layer
+	.show_images(items=[])
+		@items - if empty = all, if not - items[i] - can be layer index or layer's label name
+	.show_image(index)
+Examples:
+#1
+
+'''
+class imagej_tiff:
+	# imagej stores labels lengths in this tag
+	__TIFF_TAG_LABELS_LENGTHS = 50838
+	# imagej stores labels contents in this tag
+	__TIFF_TAG_LABELS_STRINGS = 50839
+	# init
+	def __init__(self,filename, layers = None, tile_list = None):
+		# file name
+		self.fname = filename
+		tif = Image.open(filename)
+		# total number of layers in tiff
+		self.nimages = tif.n_frames
+		# labels array
+		self.labels = []
+		# infos will contain xml data Elphel stores in some of tiff files
+		self.infos = []
+		# dictionary from decoded infos[0] xml data
+		self.props = {}
+	
+		# bits per sample, type int
+		self.bpp = tif.tag[258][0]
+
+		# Extract and parse header information
+		self.__split_labels(tif.n_frames,tif.tag)
+		self.__parse_info()
+		try:
+			self.nan_bug = self.props['VERSION']== '1.0' # data between min and max is mapped to 0..254 instead of  1.255
+		except:
+			self.nan_bug = False # other files, not ML ones
+		# image layers stacked along depth - (think RGB)
+		self.image = []
+	
+		if layers is None:
+		# fill self.image
+			for i in range(self.nimages):
+				tif.seek(i)
+				a = np.array(tif)
+				a = np.reshape(a,(a.shape[0],a.shape[1],1))
+				
+				#a = a[:,:,np.newaxis]
+				# exclude layer named 'other'
+				if self.bpp==8:
+					_min = self.data_min
+					_max = self.data_max
+					_MIN = 1
+					_MAX = 255
+					if (self.nan_bug):
+						_MIN = 0
+						_MAX = 254
+					else:
+						if self.labels[i]!='other':
+							a[a==0]=np.nan
+					a = a.astype(float)
+					if self.labels[i]!='other':
+						a = (_max-_min)*(a-_MIN)/(_MAX-_MIN)+_min
+				if i==0:
+					self.image = a
+					# stack along depth (think of RGB channels)
+				else:
+					self.image = np.append(self.image,a,axis=2)
+		else:
+			if tile_list is None:
+				indx = 0
+				for layer in layers:
+					tif.seek(self.labels.index(layer)) 
+					a = np.array(tif)
+					if not indx:
+						self.image = np.empty((a.shape[0],a.shape[1],len(layers)),a.dtype)
+					self.image[...,indx] = a
+					indx += 1
+			else:
+				other_label = "other"
+	#            print(tile_list)
+				num_tiles =  len(tile_list)
+				num_layers = len(layers)
+				tiles_corr = np.empty((num_tiles,num_layers,self.tileH*self.tileW),dtype=float)
+	#            tiles_other=np.empty((num_tiles,3),dtype=float)
+				tiles_other=self.gettilesvalues(
+						 tif = tif,
+						 tile_list=tile_list,
+						 label=other_label)
+				for nl,label in enumerate(layers):
+					tif.seek(self.labels.index(label))
+					layer = np.array(tif) # 8 or 32 bits
+					tilesX = layer.shape[1]//self.tileW
+					for nt,tl in enumerate(tile_list):
+						ty = tl // tilesX
+						tx = tl % tilesX
+	#                    tiles_corr[nt,nl] = np.ravel(layer[self.tileH*ty:self.tileH*(ty+1),self.tileW*tx:self.tileW*(tx+1)])
+						a = np.ravel(layer[self.tileH*ty:self.tileH*(ty+1),self.tileW*tx:self.tileW*(tx+1)])
+						#convert from int8
+						if self.bpp==8:
+							a = a.astype(float)
+							if np.isnan(tiles_other[nt][0]):
+								# print("Skipping NaN tile ",tl)
+								a[...] = np.nan
+							else:
+								_min = self.data_min
+								_max = self.data_max
+								_MIN = 1
+								_MAX = 255
+								if (self.nan_bug):
+									_MIN = 0
+									_MAX = 254
+								else:    
+									a[a==0] = np.nan
+								a = (_max-_min)*(a-_MIN)/(_MAX-_MIN)+_min
+						tiles_corr[nt,nl] = a    
+						pass
+					pass
+				self.corr2d =           tiles_corr
+				self.target_disparity = tiles_other[...,0]
+				self.gt_ds =            tiles_other[...,1:3]
+				pass
+		# init done, close the image
+		if (self.props['VERSION']== 2.0):
+#            self.tileH = self.image.shape[0]//self.props['tileStepY']
+#            self.tileW = self.image.shape[1]//self.props['tileStepX']
+			self.tileH = self.props['tileStepY']
+			self.tileW = self.props['tileStepX']
+			pass
+		
+		tif.close()
+	#   label == tiff layer name
+	def getvalues(self,label=""):
+		l = self.getstack([label],shape_as_tiles=True)
+		res = np.empty((l.shape[0],l.shape[1],3))
+	
+		for i in range(res.shape[0]):
+			for j in range(res.shape[1]):
+				# 9x9 -> 81x1
+				m = np.ravel(l[i,j])
+				if self.bpp==32:
+					res[i,j,0] = m[0]
+					res[i,j,1] = m[2]
+					res[i,j,2] = m[4]
+				elif self.bpp==8:
+					res[i,j,0] = ((m[0]-128)*256+m[1])/128
+					res[i,j,1] = ((m[2]-128)*256+m[3])/128
+					res[i,j,2] = (m[4]*256+m[5])/65536.0
+				else:
+					res[i,j,0] = np.nan
+					res[i,j,1] = np.nan
+					res[i,j,2] = np.nan
+	
+		# NaNize
+		a = res[:,:,0]
+		a[a==-256] = np.nan
+		b = res[:,:,1]
+		b[b==-256] = np.nan
+		c = res[:,:,2]
+		c[c==0] = np.nan
+		return res
+
+	# 3 values per tile: target disparity, GT disparity, GT confidence
+	def gettilesvalues(self,
+					 tif,
+					 tile_list,
+					 label=""):
+		res = np.empty((len(tile_list),3),dtype=float)
+		tif.seek(self.labels.index(label))
+		layer = np.array(tif) # 8 or 32 bits
+		tilesX = layer.shape[1]//self.tileW
+		for i,tl in enumerate(tile_list):
+			ty = tl // tilesX
+			tx = tl % tilesX
+			m = np.ravel(layer[self.tileH*ty:self.tileH*(ty+1),self.tileW*tx:self.tileW*(tx+1)])
+			if self.bpp==32:
+				res[i,0] = m[0]
+				res[i,1] = m[2]
+				res[i,2] = m[4]
+			elif self.bpp==8:
+				res[i,0] = ((m[0]-128)*256+m[1])/128
+				res[i,1] = ((m[2]-128)*256+m[3])/128
+				res[i,2] = (m[4]*256+m[5])/65536.0
+			else:
+				res[i,0] = np.nan
+				res[i,1] = np.nan
+				res[i,2] = np.nan
+		# NaNize
+		a = res[...,0]
+		a[a==-256] = np.nan
+		b = res[...,1]
+		b[b==-256] = np.nan
+		c = res[...,2]
+		c[c==0] = np.nan
+		
+		return res
+
+	# get ordered stack of images by provided items
+	# by index or label name. Divides into [self.tileH][self.tileW] tiles
+	def getstack(self,items=[],shape_as_tiles=False):
+		a = ()
+		if len(items)==0:
+			b = self.image
+		else:
+			for i in items:
+				if type(i)==int:
+					a += (self.image[:,:,i],)
+				elif type(i)==str:
+					j = self.labels.index(i)
+					a += (self.image[:,:,j],)
+			# stack along depth
+			b = np.stack(a,axis=2)
+		if shape_as_tiles:
+			b = get_tile_images(b,self.tileW,self.tileH)
+	
+		return b
+
+	# Trim stack to correlation data (e.g 15x15 instead of 16x16)
+	def trimStack (self, stack, radius = 0):
+		if (radius == 0):
+			radius=self.props['corrRadius']
+		corr_side = 2*radius+1
+		return stack[:,:,:,:corr_side,:corr_side]
+	# get np.array of a channel
+	# * do not handle out of bounds
+	def channel(self,index):
+		return self.image[:,:,index]
+
+
+	# extract 2D corr data
+	def get2DCorr(self):
+		im_2DCorr = crop_border(self.image, self.tileH)
+		#print('im_2DCorr shape', im_2DCorr.shape)
+		im_2DCorr = delete_tile_splits(im_2DCorr, self.tileH)
+		#print('im_2DCorr shape', im_2DCorr.shape)
+		# Delete last 2 layers
+		# im_2DCorr = np.delete(im_2DCorr,im_2DCorr.shape[2]-1,axis=2)
+		im_2DCorr = np.delete(im_2DCorr,[120,121],axis=2)
+		return im_2DCorr
+
+
+	# Retrieve specific meta image defined by index & layer
+	def getMeta(self, index, layer = 0):
+		_, meta = self.getCorrsMeta(items=[layer])
+		meta = np.squeeze(meta)
+		return meta[:,:,index]
+
+
+	# Retrieve correlation information (tiles) and meta-data
+	def getCorrsMeta(self,items=[]):
+		stack0 =  self.getstack(items,shape_as_tiles=True)
+		stack = np.moveaxis(stack0, 4, 0) # slices - first index
+		
+		radius=self.props['corrRadius']
+		num_meta=self.props['numMeta']
+		corr_side = 2*radius+1
+		corr_tiles = stack[:,:,:,:corr_side,:corr_side]
+		
+		meta = stack[:,:,:,-1,:num_meta]
+		return corr_tiles, meta/self.props['tileMetaScale'] 
+
+	# display images by index or label
+	def show_images(self,items=[]):
+
+		# show listed only
+		if len(items)>0:
+			for i in items:
+				if type(i)==int:
+					self.show_image(i)
+				elif type(i)==str:
+					j = self.labels.index(i)
+					self.show_image(j)
+		# show all
+		else:
+			for i in range(self.nimages):
+				self.show_image(i)
+
+
+	# display single image
+	def show_image(self,index):
+	# display using matplotlib
+
+		t = self.image[:,:,index]
+		mytitle = "("+str(index+1)+" of "+str(self.nimages)+") "+self.labels[index]
+		fig = plt.figure()
+		fig.canvas.set_window_title(self.fname+": "+mytitle)
+		fig.suptitle(mytitle)
+		#plt.imshow(t,cmap=plt.get_cmap('gray'))
+		plt.imshow(t)
+		plt.colorbar()
+	
+		# display using Pillow - need to scale
+	
+		# remove NaNs - no need
+		#t[np.isnan(t)]=np.nanmin(t)
+		# scale to [min/max*255:255] range
+		#t = (1-(t-np.nanmax(t))/(t-np.nanmin(t)))*255
+		#tmp_im = Image.fromarray(t)
+		#tmp_im.show()
+
+	# puts etrees in infos
+	def __parse_info(self):
+
+		infos = []
+		for info in self.infos:
+			infos.append(ET.fromstring(info))
+		self.infos = infos
+	
+		# specifics
+		# properties dictionary
+		pd = {}
+	
+		if infos:
+			for child in infos[0]:
+				#print(child.tag+"::::::"+child.text)
+				pd[child.tag] = child.text
+		
+			self.props = pd
+			file_version = float(self.props['VERSION'])
+			if (file_version < 2.0):
+				# tiles are squares (older version
+				self.tileW = int(self.props['tileWidth'])
+				self.tileH = int(self.props['tileWidth'])
+				self.data_min = float(self.props['data_min'])
+				self.data_max = float(self.props['data_max'])
+			else:
+				floats=['dispOffsetLow','tileMetaScale','disparity_low','dispOffset',
+				'fatZero','disparity_pwr','VERSION','dispOffsetHigh',
+				'disparity_high']
+				ints = ['metaGTConfidence','tileMetaSlice','indexReference','metaLastDiff',
+				 'metaGTDisparity','metaFracValid', 'numScenes','metaTargetDisparity',
+				 'disparity_steps','tileStepX',    'tileStepY', "corrRadius","numMeta"]
+				bools=['randomize_offsets']
+				for key in pd:
+					val = pd[key]
+					if key in bools:
+						if (val == '1') or (val == 'true') or (val == 'True'):
+							pd[key] = 1
+						else:
+							pd[key] = 0  
+						pass
+					elif key in ints:
+						pd[key] = int(pd[key])
+					elif  key in floats:   
+						pd[key] = float(pd[key])
+				try:       
+					pd['corrRadius'] = pd['corrRadius'] # not yet exists
+				except:
+					pd['corrRadius'] = 7              
+				try:       
+					pd['numMeta'] = pd['numMeta'] # not yet exists
+				except:
+					pd['numMeta'] = 6              
+						
+			pass        
+
+	# makes arrays of labels (strings) and unparsed xml infos
+	def __split_labels(self,n,tag):
+		# list
+		tag_lens = tag[self.__TIFF_TAG_LABELS_LENGTHS]
+		# string
+		tag_labels = tag[self.__TIFF_TAG_LABELS_STRINGS].decode()
+		# remove 1st element: it's something like IJIJlabl..
+		tag_labels = tag_labels[tag_lens[0]:]
+		tag_lens = tag_lens[1:]
+		# the last ones are images labels
+		# normally the difference is expected to be 0 or 1
+		skip = len(tag_lens) - n
+	
+		self.labels = []
+		self.infos = []
+		for l in tag_lens:
+			string = tag_labels[0:l].replace('\x00','')
+			if skip==0:
+				self.labels.append(string)
+			else:
+				self.infos.append(string)
+				skip -= 1
+			tag_labels = tag_labels[l:]
+
+
+def arg_parser():
+	parser = argparse.ArgumentParser(description='Extract images from specific multi-layer tiff file')
+	required = parser.add_argument_group('Required')
+	required.add_argument('-i', '--input', type=str, required=True,
+						  help='input TIFF file (specific multi-layer)')
+	options = parser.add_argument_group('Options')
+	options.add_argument('--crop_border', action="store_true",
+						  help='Crop image border (removing NaN border)')
+	options.add_argument('--corr', type=str,
+						  help='Save 2D Corr TIFF file (multi-layer)')
+	options.add_argument('--targetdisp', type=str,
+						  help='Save Target Disparity TIFF file (single layer)')
+	options.add_argument('--groundtruth', type=str,
+						  help='Save Ground Truth TIFF file (single layer)')
+	options.add_argument('--confidence', type=str, 
+						  help='Save ground truth confidence image (single layer')
+	options.add_argument('--disp_lma', type=str, 
+						  help='Save LMA disparity image (single layer')
+	options.add_argument('--frac_valid', type=str, 
+						  help='Save fraction valid image (single layer')
+	options.add_argument('--last_diff', type=str, 
+						  help='Save image (single layer')
+	options.add_argument('--disp_bg', type=str, 
+						  help='Save image (single layer')
+	options.add_argument('--confidence_bg', type=str, 
+						  help='Save image (single layer')
+	options.add_argument('--disp_lma_bg', type=str, 
+						  help='Save image (single layer')
+	options.add_argument('--last_diff_bg', type=str, 
+						  help='Save image (single layer')
+	options.add_argument('--disp_fg', type=str, 
+						  help='Save image (single layer')
+	options.add_argument('--disp_bg_all', type=str, 
+						  help='Save image (single layer')
+	options.add_argument('--blue_sky', type=str, 
+						  help='Save image (single layer')
+	options.add_argument('--verbose', action="store_true",
+						  help='verbose mode')
+	return parser
+
+
+#MAIN
+def main(args=None):
+	args = arg_parser().parse_args(args)
+
+	Time1 = time.time()
+	#print(bcolors.BOLDWHITE+"Time1: "+str(Time1)+bcolors.ENDC)
+
+	print('Reading Tiff image...')
+	ijt = imagej_tiff(args.input)
+
+	Time2 = time.time()
+	#print(bcolors.BOLDWHITE+"Time2: "+str(Time2)+bcolors.ENDC)
+	TimeDiff = Time2 - Time1
+	print(bcolors.BOLDWHITE+"Reading Time: "+str(TimeDiff)+bcolors.ENDC)
+
+	# Image header info
+	rough_string = ET.tostring(ijt.infos[0], "utf-8")
+	reparsed = minidom.parseString(rough_string)
+	print('\n Image Header info:')
+	print(reparsed.toprettyxml(indent="\t"))
+
+	print("\n TIFF stack labels: "+str(ijt.labels))
+	print("TIFF stack labels length: "+str(len(ijt.labels)))
+	#print(ijt.infos)
+
+
+	# Image shape information and tiles information
+	print("\nImage shape: ",ijt.image.shape)
+	print("Image size: ",ijt.image.size)
+	print("Image strides: ",ijt.image.strides)
+
+	# needed properties:
+	print("Tiles shape: "+str(ijt.tileW)+"x"+str(ijt.tileH))
+
+	try:
+		print("Data min: "+str(ijt.data_min))
+		print("Data max: "+str(ijt.data_max))
+	except:
+		print("Data min/max are not provided")
+
+
+#    tiles,tiles_meta = ijt.getCorrsMeta(['0-1','1-2','2-3','3-4'])
+	# nrows, ncols, height, width, depth
+	tiles,tiles_meta = ijt.getCorrsMeta([])
+	print("Corr stack shape: "+str(tiles.shape))
+	print("Meta stack shape: "+str(tiles_meta.shape))
+
+
+	# Reorder numpy array
+	# img_reordered = np.transpose(ijt.image, axes=[2,0,1])
+	# print('img_reordered shape',img_reordered.shape)
+	# imageio.mimwrite('test_imageio.tiff',img_reordered)
+	# #tifffile.imwrite('test_tifffile.tiff',img_reordered)
+
+
+	# Extract and crop 2d correlation data
+	if (args.corr is not None):
+		Corr_FileName = args.corr
+		im_2DCorr = ijt.get2DCorr()
+		im_2DCorr = np.moveaxis(im_2DCorr, 2, 0)
+		if args.verbose:
+			print("\nImage shape: ",ijt.image.shape)
+			print('im_2DCorr shape', im_2DCorr.shape)
+		print('\nSaving 2DCorr image... ')
+		if not os.path.exists(os.path.dirname(Corr_FileName)):
+			os.makedirs(os.path.dirname(Corr_FileName))
+		imageio.mimwrite(Corr_FileName, im_2DCorr)
+
+	# Extract and crop target disparity
+	if (args.targetdisp is not None):
+		#TargetDisparity = tiles_meta[0,:,:,0]
+		save_metadata_layer(ijt, args.targetdisp, 0, args.crop_border, args.verbose)
+		print('Saving targetdisp layer...')
+	
+	# Extract and crop Ground Truth disparity
+	if (args.groundtruth is not None):
+		#GroundTruth = tiles_meta[0,:,:,1]
+		save_metadata_layer(ijt, args.groundtruth, 1, args.crop_border, args.verbose)
+		print('Saving groundtruth layer...')
+		
+	# Extract and crop ground truth confidence
+	if (args.confidence is not None):
+		#Confidence = tiles_meta[0,:,:,2]
+		save_metadata_layer(ijt, args.confidence, 2, args.crop_border, args.verbose)
+		print('Saving confidence layer...')
+	
+	# Extract and crop LMA disparity image
+	if (args.disp_lma is not None):
+		#DispLMA = tiles_meta[0,:,:,3]
+		save_metadata_layer(ijt, args.disp_lma, 3, args.crop_border, args.verbose)
+		print('Saving disp_lma layer...')
+	
+	# Extract and crop frac_valid
+	if (args.frac_valid is not None):
+		#FracValid = tiles_meta[0,:,:,4]
+		save_metadata_layer(ijt, args.frac_valid, 4, args.crop_border, args.verbose)
+		print('Saving frac_valid layer...')
+	
+	# Extract and crop last_diff
+	if (args.last_diff is not None):
+		#LastDiff = tiles_meta[0,:,:,5]
+		save_metadata_layer(ijt, args.last_diff, 5, args.crop_border, args.verbose)
+		print('Saving last_diff layer...')
+	
+	# Extract and crop disp_bg
+	if (args.disp_bg is not None):
+		#DispBg = tiles_meta[0,:,:,6]
+		save_metadata_layer(ijt, args.disp_bg, 6, args.crop_border, args.verbose)
+		print('Saving disp_bg layer...')
+	
+	# Extract and crop confidence_bg
+	if (args.confidence_bg is not None):
+		#ConfidenceBg = tiles_meta[0,:,:,7]
+		save_metadata_layer(ijt, args.confidence_bg, 7, args.crop_border, args.verbose)
+		print('Saving confidence_bg layer...')
+	
+	# Extract and crop disp_lma_bg
+	if (args.disp_lma_bg is not None):
+		#DispLMABg = tiles_meta[0,:,:,8]
+		save_metadata_layer(ijt, args.disp_lma_bg, 8, args.crop_border, args.verbose)
+		print('Saving disp_lma_bg layer...')
+	
+	# Extract and crop last_diff_bg
+	if (args.last_diff_bg is not None):
+		#LastDiffBg = tiles_meta[0,:,:,9]
+		save_metadata_layer(ijt, args.last_diff_bg, 9, args.crop_border, args.verbose)
+		print('Saving last_diff_bg layer...')
+	
+	# Extract and crop disp_fg
+	if (args.disp_fg is not None):
+		#DispFg = tiles_meta[0,:,:,10]
+		save_metadata_layer(ijt, args.disp_fg, 10, args.crop_border, args.verbose)
+		print('Saving disp_fg layer...')
+	
+	# Extract and crop disp_bg_all
+	if (args.disp_bg_all is not None):
+		#DispBgAll = tiles_meta[0,:,:,11]
+		save_metadata_layer(ijt, args.disp_bg_all, 11, args.crop_border, args.verbose)
+		print('Saving disp_bg_all layer...')
+	
+	# Extract and crop blue_sky
+	if (args.blue_sky is not None):
+		#BlueSky = tiles_meta[0,:,:,12]
+		save_metadata_layer(ijt, args.blue_sky, 12, args.crop_border, args.verbose)
+		print('Saving blue_sky layer...')
+	
+	Time3 = time.time()
+	#print(bcolors.BOLDWHITE+"Time2: "+str(Time2)+bcolors.ENDC)
+	TimeDiff = Time3 - Time1
+	print(bcolors.BOLDWHITE+"Execution Time: "+str(TimeDiff)+bcolors.ENDC)
+
+
+if __name__ == '__main__':
+	sys.exit(main(sys.argv[1:]))

Preprocessing_CombinedImages.py

+#!/usr/bin/env python3
+
+
+import imageio 
+import numpy as np
+
+import sys
+import argparse
+import time 
+
+def arg_parser():
+	parser = argparse.ArgumentParser(description='Combine images to generate 3D volume')
+	required = parser.add_argument_group('Required')
+	required.add_argument('--corr', type=str, required=True,
+						  help='2D Corr TIFF file (multi-layer)')
+	required.add_argument('--targetdisp', type=str, required=True,
+						  help='Target Disparity TIFF file (single layer)')
+	required.add_argument('--groundtruth', type=str, required=True,
+						  help='Ground Truth TIFF file (single layer)')
+	required.add_argument('--confidence', type=str, required=True,
+						  help='Confidence TIFF file (single layer)')
+	required.add_argument('--disp_lma', type=str, required=True,
+						  help='LMA Disparity TIFF file (single layer)')
+	required.add_argument('--output', type=str, required=True,
+						  help='Combined TIFF file (multi-layer)')
+	options = parser.add_argument_group('Options')
+	options.add_argument('--repeat', type=int, default=15,
+						  help='Repeat tile number (default 15)')
+	return parser
+
+
+#MAIN
+def main(args=None):
+	args = arg_parser().parse_args(args)
+
+	Time1 = time.time()
+	#print(bcolors.BOLDWHITE+"Time1: "+str(Time1)+bcolors.ENDC)
+
+	img2DCorr_name = args.corr
+	imgTargetDisp_name = args.targetdisp
+	imgGT_name = args.groundtruth
+	imgConfidence_name = args.confidence
+	imgDispLMA_name = args.disp_lma
+	output_name = args.output
+	repeatNb = args.repeat
+	
+	# Read 2dcorr image - all layers
+	print("Reading 2dcorr file...")
+	img2DCorr = imageio.mimread(img2DCorr_name,memtest=False)
+	img2DCorr = np.array(img2DCorr)
+	# print('\nimg2DCorr type: ', img2DCorr.dtype)
+	# print('img2DCorr shape: ', img2DCorr.shape)
+
+	# Read TargetDisp and GroundTruth images
+	print("Reading TargetDisp file...")
+	imgTargetDisp = imageio.imread(imgTargetDisp_name)
+	# print('imgTargetDisp type: ', imgTargetDisp.dtype)
+	# print('imgTargetDisp shape: ', imgTargetDisp.shape)
+
+	print("Reading GroundTruth file...")
+	imgGT = imageio.imread(imgGT_name)
+	# print('imgGT type: ', imgGT.dtype)
+	# print('imgGT shape: ', imgGT.shape)
+
+	print("Reading Confidence file...")
+	imgConfidence = imageio.imread(imgConfidence_name)
+	# print('imgConfidence type: ', imgConfidence.dtype)
+	# print('imgConfidence shape: ', imgConfidence.shape)
+
+	print("Reading DispLMA file...")
+	imgDispLMA = imageio.imread(imgDispLMA_name)
+	# print('imgDispLMA type: ', imgDispLMA.dtype)
+	# print('imgDispLMA shape: ', imgDispLMA.shape)
+
+	# - - - - - - - - - - -
+	print("Generating combined image...")
+	# Resample TargetDisp and imgTG (repeating values, to match img2DCor size)
+	imgTargetDisp_repeat0 = np.repeat(imgTargetDisp, repeatNb, axis=0)
+	imgTargetDisp_repeat = np.repeat(imgTargetDisp_repeat0, repeatNb, axis=1)
+	imgTargetDisp_repeat = np.expand_dims(imgTargetDisp_repeat, axis=0)
+
+	imgGT_repeat0 = np.repeat(imgGT, repeatNb, axis=0)
+	imgGT_repeat = np.repeat(imgGT_repeat0, repeatNb, axis=1)
+	imgGT_repeat = np.expand_dims(imgGT_repeat, axis=0)
+	# print('imgTargetDisp_repeat shape: ', imgTargetDisp_repeat.shape)
+	# print('imgGT_repeat shape: ', imgGT_repeat.shape)
+
+	imgConfidence_repeat0 = np.repeat(imgConfidence, repeatNb, axis=0)
+	imgConfidence_repeat = np.repeat(imgConfidence_repeat0, repeatNb, axis=1)
+	imgConfidence_repeat = np.expand_dims(imgConfidence_repeat, axis=0)
+	print('\t imgConfidence_repeat shape: ', imgConfidence_repeat.shape)
+
+	imgDispLMA_repeat0 = np.repeat(imgDispLMA, repeatNb, axis=0)
+	imgDispLMA_repeat = np.repeat(imgDispLMA_repeat0, repeatNb, axis=1)
+	imgDispLMA_repeat = np.expand_dims(imgDispLMA_repeat, axis=0)
+	print('\t imgDispLMA_repeat shape: ', imgDispLMA_repeat.shape)
+
+	# Stack layers to img2dCorr, to generate one 3D volume
+	imgAll = np.concatenate((img2DCorr,imgTargetDisp_repeat), axis = 0)
+	imgAll = np.concatenate((imgAll,imgGT_repeat), axis = 0)
+	imgAll = np.concatenate((imgAll,imgConfidence_repeat), axis = 0)
+	imgAll = np.concatenate((imgAll,imgDispLMA_repeat), axis = 0)
+	
+	imageio.mimwrite(output_name,imgAll)
+
+	Time2 = time.time()
+	#print(bcolors.BOLDWHITE+"Time2: "+str(Time2)+bcolors.ENDC)
+	TimeDiff = Time2 - Time1
+	print("Reading Time: "+str(TimeDiff))
+
+if __name__ == '__main__':
+	sys.exit(main(sys.argv[1:]))
+

QualityControl_Confidence.py

+#!/usr/bin/env python3
+
+
+import imageio 
+import numpy as np
+import pandas as pd
+
+import sys
+import argparse
+import time 
+import os 
+
+def arg_parser():
+	parser = argparse.ArgumentParser(description='Quality control')
+	required = parser.add_argument_group('Required')
+	required.add_argument('--confidence', type=str, required=True,
+						  help='2D confidence TIFF file (single layer)')
+	required.add_argument('--output', type=str, required=True,
+						  help='CSV file')
+	return parser
+
+
+#MAIN
+def main(args=None):
+	args = arg_parser().parse_args(args)
+
+	Time1 = time.time()
+	#print(bcolors.BOLDWHITE+"Time1: "+str(Time1)+bcolors.ENDC)
+
+	imgConfidence_name = args.confidence
+	imgConfidence_basename = os.path.basename(imgConfidence_name)
+	output_name = args.output
+
+	# Read confidence image - all layers
+	#print("Reading confidence file...")
+	imgConfidence = imageio.imread(imgConfidence_name)
+	print('\nimgConfidence type: ', imgConfidence.dtype)
+	print('imgConfidence shape: ', imgConfidence.shape)
+
+	Bool_NaN = np.isnan(imgConfidence)
+	Nb_NaN = np.count_nonzero(Bool_NaN)
+	imgConfidence_NotNaN = imgConfidence[~np.isnan(imgConfidence)]
+	print('\nNb_NaN: ', Nb_NaN)
+	print('\nimgConfidence_NotNaN type: ', imgConfidence_NotNaN.dtype)
+	print('imgConfidence_NotNaN shape: ', imgConfidence_NotNaN.shape)
+
+	# Compute Image info: mean,min,max values
+	imgConfidence_mean = np.mean(imgConfidence_NotNaN)
+	imgConfidence_min = np.min(imgConfidence_NotNaN)
+	imgConfidence_max = np.max(imgConfidence_NotNaN)	
+	print('\nimgConfidence_mean : ', imgConfidence_mean)
+	print('imgConfidence_min : ', imgConfidence_min)
+	print('imgConfidence_max : ', imgConfidence_max)
+	
+	# Save information to CSV file
+	QC_data = np.array([[imgConfidence_basename,imgConfidence.shape[0],imgConfidence.shape[1], \
+		Nb_NaN, imgConfidence_mean, imgConfidence_min, imgConfidence_max]
+	])
+	Columns = ['FileName','imgConfidence_Shape0','imgConfidence_Shape1', \
+	'Nb_NaN','imgConfidence_mean','imgConfidence_min','imgConfidence_max']
+	df = pd.DataFrame(QC_data,columns=Columns)
+	print(df.head())
+	df.to_csv(output_name, index=False)
+
+	Time2 = time.time()
+	#print(bcolors.BOLDWHITE+"Time2: "+str(Time2)+bcolors.ENDC)
+	TimeDiff = Time2 - Time1
+	#print("Computing Time: "+str(TimeDiff))
+
+if __name__ == '__main__':
+	sys.exit(main(sys.argv[1:]))
+

QualityControl_Images.py

+#!/usr/bin/env python3
+
+
+import imageio 
+import numpy as np
+import pandas as pd
+
+import sys
+import argparse
+import time 
+import os 
+
+def arg_parser():
+	parser = argparse.ArgumentParser(description='Quality control')
+	required = parser.add_argument_group('Required')
+	required.add_argument('--corr', type=str, required=True,
+						  help='2D Corr TIFF file (multi-layer)')
+	required.add_argument('--targetdisp', type=str, required=True,
+						  help='Target Disparity TIFF file (single layer)')
+	required.add_argument('--groundtruth', type=str, required=True,
+						  help='Ground Truth TIFF file (single layer)')
+	required.add_argument('--output', type=str, required=True,
+						  help='CSV file')
+	return parser
+
+
+#MAIN
+def main(args=None):
+	args = arg_parser().parse_args(args)
+
+	Time1 = time.time()
+	#print(bcolors.BOLDWHITE+"Time1: "+str(Time1)+bcolors.ENDC)
+
+	img2DCorr_name = args.corr
+	img2DCorr_basename = os.path.basename(img2DCorr_name)
+	imgTargetDisp_name = args.targetdisp
+	imgGT_name = args.groundtruth
+	output_name = args.output
+
+	# Read 2dcorr image - all layers
+	#print("Reading 2dcorr file...")
+	img2DCorr = imageio.mimread(img2DCorr_name,memtest=False)
+	img2DCorr = np.array(img2DCorr)
+	# print('\nimg2DCorr type: ', img2DCorr.dtype)
+	#print('img2DCorr shape: ', img2DCorr.shape)
+
+	# Read TargetDisp and GroundTruth images
+	#print("Reading TargetDisp file...")
+	imgTargetDisp = imageio.imread(imgTargetDisp_name)
+	# print('imgTargetDisp type: ', imgTargetDisp.dtype)
+	#print('imgTargetDisp shape: ', imgTargetDisp.shape)
+
+	#print("Reading GroundTruth file...")
+	imgGT = imageio.imread(imgGT_name)
+	# print('imgGT type: ', imgGT.dtype)
+	#print('imgGT shape: ', imgGT.shape)
+
+	# Compute Image info: mean,min,max values
+	img2DCorr_mean = np.mean(img2DCorr)
+	img2DCorr_min = np.min(img2DCorr)
+	img2DCorr_max = np.max(img2DCorr)
+	
+	imgTargetDisp_mean = np.mean(imgTargetDisp)
+	imgTargetDisp_min = np.min(imgTargetDisp)
+	imgTargetDisp_max = np.max(imgTargetDisp)
+	
+	imgGT_mean = np.mean(imgGT)
+	imgGT_min = np.min(imgGT)
+	imgGT_max = np.max(imgGT)
+	# print('\n img2DCorr pixel info:')
+	# print('\t img2DCorr_mean: ', img2DCorr_mean)
+	# print('\t img2DCorr_min: ', img2DCorr_min)
+	# print('\t img2DCorr_max: ', img2DCorr_max)
+
+	# Compute image difference between TargetDisp and GroundTruth images
+	Diff = np.abs(imgTargetDisp - imgGT)
+	Diff_mean = np.mean(Diff)
+	Diff_min = np.min(Diff)
+	Diff_max = np.max(Diff)
+	# print('Image difference:')
+	# print('\t Diff_mean:',Diff_mean )
+	# print('\t Diff_min:', Diff_min)
+	# print('\t Diff_max:', Diff_max)
+	
+
+	# Save information to CSV file
+	QC_data = np.array([[img2DCorr_basename,img2DCorr.shape[0],img2DCorr.shape[1], img2DCorr.shape[2], \
+		imgTargetDisp.shape[0], imgTargetDisp.shape[1], imgGT.shape[0], imgGT.shape[1], \
+		img2DCorr_mean, img2DCorr_min, img2DCorr_max, \
+		imgTargetDisp_mean, imgTargetDisp_min, imgTargetDisp_max, \
+		imgGT_mean, imgGT_min, imgGT_max, \
+		Diff_mean, Diff_min, Diff_max]
+	])
+	Columns = ['FileName','img2DCorr_Shape0','img2DCorr_Shape1','img2DCorr_Shape2',\
+	'imgTargetDisp_Shape0','imgTargetDisp_Shape1',\
+	'imgGT_Shape0','imgGT_Shape1',
+	'img2DCorr_mean','img2DCorr_min','img2DCorr_max',\
+	'imgTargetDisp_mean','imgTargetDisp_min','imgTargetDisp_max',\
+	'imgGT_mean','imgGT_min','imgGT_max',\
+	'Diff_mean','Diff_min','Diff_max',\
+	]
+	df = pd.DataFrame(QC_data,columns=Columns)
+	print(df.head())
+	df.to_csv(output_name, index=False)
+
+	Time2 = time.time()
+	#print(bcolors.BOLDWHITE+"Time2: "+str(Time2)+bcolors.ENDC)
+	TimeDiff = Time2 - Time1
+	#print("Computing Time: "+str(TimeDiff))
+
+if __name__ == '__main__':
+	sys.exit(main(sys.argv[1:]))
+

README.txt

+This project contains various python scripts used for disparity map analysis
+
+
+## 1. Data pre-processing stage
+
+Note: Data preprocessing stage to generate input files for AI analysis via ir-tp-net
+
+### 1.1. Extracting images of interest from multi-layer 3D TIFF files
+
+Example:
+''' 
+python3 ./Preprocessing_2DPhaseCorrelation.py --crop_border --input $i --targetdisp $TargetDispImage \
+		--groundtruth $GroundTruthImage --confidence $ConfidenceImage --disp_lma $DispLMAImage --corr $CorrImage --verbose
+'''
+
+Notes:
+ - Input file is a special multi-layer TIFF file
+ - Output CorrImage is a 3D TIFF file with 120 layers
+ - Other output files are 2D images
+
+### 1.2. Generating 3D TIFF files as direct input files to neural network ir-tp-net
+
+Example:
+''' 
+python3 ./Preprocessing_CombinedImages.py --corr $CorrImage --targetdisp $TargetDispImage \
+		--groundtruth $GroundTruthImage --confidence $ConfidenceImage --disp_lma $DispLMAImage --output $CombinedImage
+''' 
+Notes: 
+ - Output file is a 3D image with 124 layers
+  - scaled to the size of the correlation image
+
+## 2. AI analysis stage
+
+Please see "ir-tp-net" project for neural network training and testing to generate predicted disparity map
+
+
+## 3. Data post-processing stage
+
+Note: Data postprocessing stage to analyze outputs from AI analysis
+
+###  3.1. Density analysis
+
+Examples:
+''' 
+python3 ./Compute_Density.py --pred $i --groundtruth $GroundTruthImage --adjtilesdim 1 --output $Density_CSVFile --inclusionmask $InclusionMask_File --exclusionmask $ExclusionMask_File --threshold 2.0 --verbose
+''' 
+
+Note: input file is the predicted disparity map
+
+###  3.2. RMSE analysis
+
+Example:
+''' 
+python3 ./Compute_RMSE_WithThreshold.py --pred $i --groundtruth $GroundTruthImage --adjtilesdim $AdjTilesDim --threshold $RMSE_Threshold --output $RMSE_CSVFile
+			
+python3 ./Compute_RMSE_WithFiltering.py --pred $i --groundtruth $GroundTruthImage --confidence $ConfidenceImage --disp_lma $DispLMAImage --adjtilesdim 1 --threshold $Threshold --output $RMSE_CSVFile --output_mask $RMSEFiltering_MaskFile
+''' 
+
+Note: input file is the predicted disparity map
+
+
+## 4. Data Quality control
+
+Quality control stage generating multi-layer 3D tiff file
+
+''' 
+python3 ./Compute_Inference_QCImage.py --pred $i --groundtruth $GroundTruthImage --targetdisp $TargetDispImage --mask $DataFilteringMask --threshold 2.0 --output $InferenceQC_File --verbose
+''' 
+
+Notes: 
+ - input file is the predicted disparity map
+ - output file includes predicted disparity map, ground truth map, target disparity map and mask map

RandomNoiseLevel_ImageCreation.py

+#!/usr/bin/env python3
+
+
+import imageio 
+import numpy as np
+
+import sys
+import argparse
+import time 
+
+def arg_parser():
+	parser = argparse.ArgumentParser(description='Generate combined image with random noise level')
+	required = parser.add_argument_group('Required')
+	required.add_argument('--template', type=str, required=True,
+						  help='Combined TIFF file (multi-layer)')
+	required.add_argument('--offsetTemplate', type=float, required=True,
+						  help='Template noise offset level')
+	required.add_argument('--noiseLevel', type=float, required=True,
+						  help='Noise offset level')
+	required.add_argument('--output', type=str, required=True,
+						  help='Output combined TIFF file (multi-layer)')
+	options = parser.add_argument_group('Options')
+	options.add_argument('--tileSize', type=int, default=15,
+						  help='Tile Size')
+	options.add_argument('--outputTargetDisp', type=str,
+						  help='Output TargetDisp (2D image)')
+	options.add_argument('--outputGroundTruth', type=str,
+						  help='Output GroundTruth (2D image)')
+	return parser
+
+
+# Steps
+# Option to define offset noise level
+# Assess offset images based on noise level
+# Load all images as list
+# for each tile, pick among available images
+# Save image
+
+#MAIN
+def main(args=None):
+	args = arg_parser().parse_args(args)
+
+	Time1 = time.time()
+	#print(bcolors.BOLDWHITE+"Time1: "+str(Time1)+bcolors.ENDC)
+
+	template_name = args.template
+	offset_template = args.offsetTemplate
+	noise_level = args.noiseLevel
+	output_name = args.output
+	tile_size = args.tileSize
+	output_TargetDisp = args.outputTargetDisp
+	output_GroundTruth = args.outputGroundTruth
+
+	# Offsets
+	List_AvailableOffsets = np.array([-5.000, -4.003, -3.116, -2.341, -1.676, -1.122, -0.679, -0.346, -0.125, -0.014, \
+		0.014, 0.125, 0.346, 0.679, 1.122, 1.676, 2.341, 3.116, 4.003, 5.000])
+	#print('List_AvailableOffsets',List_AvailableOffsets)
+
+	# Filter offsets based on offsetLevel
+	List_Offsets = List_AvailableOffsets[np.logical_and(List_AvailableOffsets >= -noise_level,List_AvailableOffsets <= noise_level)]
+	print('\n List_Offsets',List_Offsets)
+
+	# Nb images
+	nb_images = len(List_Offsets)
+	print('nb_images',nb_images)
+
+	# Read template combined image
+	print("\n Reading template file...")
+	imgTemplate_layerList = imageio.mimread(template_name,memtest=False)
+	imgTemplate = np.stack(imgTemplate_layerList, axis=0)
+	print('\t imgTemplate type: ', imgTemplate.dtype)
+	print('\t imgTemplate shape: ', imgTemplate.shape)
+
+	print("\n Creating image stack...")
+	image_stack = np.zeros((nb_images, imgTemplate.shape[0], imgTemplate.shape[1], imgTemplate.shape[2]), dtype = imgTemplate.dtype)
+	# print('image_stack type: ', image_stack.dtype)
+	# print('image_stack shape: ', image_stack.shape)
+
+	# List_ImageNames
+	List_ImageNames = []
+	for i, offset in enumerate(List_Offsets):
+		print('\t Reading Image {} with offset {:.3f} ...'.format(i,offset))
+		offset_string = "{:.3f}".format(offset)
+		FileName_current = template_name.replace(str(offset_template),offset_string)
+		List_ImageNames.append(FileName_current)
+		img_current_layerList = imageio.mimread(FileName_current,memtest=False)
+		img_current = np.stack(img_current_layerList, axis=0)
+		image_stack[i,...] = img_current
+	#print('List_ImageNames',List_ImageNames)
+
+
+	# Generating image stack
+	# image_stack = np.stack(List_Images, axis=0)
+	print('\t image_stack type: ', image_stack.dtype)
+	print('\t image_stack shape: ', image_stack.shape)
+
+	# Output image
+	print("\n Generating output image by random input patch selection...")
+	imgOutput = np.zeros(imgTemplate.shape, dtype = imgTemplate.dtype)
+	print('\t imgOutput type: ', imgOutput.dtype)
+	print('\t imgOutput shape: ', imgOutput.shape)
+
+	i = j = 0
+	for i in range(0,imgTemplate.shape[1], tile_size):
+		for j in range(0,imgTemplate.shape[2], tile_size):
+			random_selection = np.random.randint(nb_images)
+			#print('i j rnd: {} {} {}'.format(i,j,random_selection))
+			#imgOutput[:,i:i+tile_size,j:j+tile_size] = tile_size * patchNb
+			imgOutput[:,i:i+tile_size,j:j+tile_size] = image_stack[random_selection,:,i:i+tile_size,j:j+tile_size]
+
+	imgOutput_TargetDisp = imgOutput[-2,:,:]
+	imgOutput_TargetDisp = imgOutput_TargetDisp[::tile_size,::tile_size]
+	print('\t imgOutput_TargetDisp shape: ', imgOutput_TargetDisp.shape)
+
+	imgOutput_GroundTruth = imgOutput[-1,:,:]
+	imgOutput_GroundTruth = imgOutput_GroundTruth[::tile_size,::tile_size]
+	print('\t imgOutput_GroundTruth shape: ', imgOutput_GroundTruth.shape)
+
+	# Write output combined image (multi-layer)
+	print("\n Writing output file - combined image...")
+	imageio.mimwrite(output_name,imgOutput)
+
+	# Write output TargetDisp image (2D)
+	if (output_TargetDisp is not None):
+		print("\nWriting output file - TargetDisp image...")
+		imageio.imwrite(output_TargetDisp,imgOutput_TargetDisp)
+
+	# Write output GroundTruth image (2D)
+	if (output_GroundTruth is not None):
+		print("\nWriting output file - GroundTruth image...")
+		imageio.imwrite(output_GroundTruth,imgOutput_GroundTruth)
+
+	Time2 = time.time()
+	#print(bcolors.BOLDWHITE+"Time2: "+str(Time2)+bcolors.ENDC)
+	TimeDiff = Time2 - Time1
+	print("Reading Time: "+str(TimeDiff))
+
+if __name__ == '__main__':
+	sys.exit(main(sys.argv[1:]))
+