diff --git a/src/argpopt.h b/src/argpopt.h
index 1a40e30..7bee75a 100644
--- a/src/argpopt.h
+++ b/src/argpopt.h
@@ -40,7 +40,7 @@ struct Config
 const char *argp_program_version = "0.2";
 const char *argp_program_bug_address = "<carl@uvos.xyz>";
 static char doc[] = "Program to determine the lubricant thikness on a curved surface.\n\
-Possible operations: apply create curve mkcurve mkboard";
+Possible operations: apply create curve mkcurve mkboard show";
 static char args_doc[] = "[OPERATION] IMAGE1 IMAGE2 ...";
 
 static struct argp_option options[] = 
diff --git a/src/main.cpp b/src/main.cpp
index d1395e8..23d6265 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -42,6 +42,7 @@ enum {
 	APPLY_MAP,
 	APPLY_CURVE,
 	CREATE_CURVE,
+	SHOW_IMAGE,
 	EXIT
 };
 
@@ -59,6 +60,8 @@ int selectOperation(char** opt)
 		return CREATE_CURVE;
 	else if(strcmp(opt[0], "mkboard" ) == 0)
 		return CREATE_CHARUCO;
+	else if(strcmp(opt[0], "show" ) == 0)
+		return SHOW_IMAGE;
 	else if(strcmp(opt[0], "exit" ) == 0)
 		return EXIT;
 	return -1;
@@ -115,13 +118,19 @@ std::vector<cv::Mat> loadImages(char** fileNames)
 		cv::Mat tmpImage;
 		const std::string str(fileNames[i]);
 		if(str.find(".mat") !=  std::string::npos)
-			openImageYml(fileNames[i]);
-		else	
-		 tmpImage = openImageImg(fileNames[i]);
+		{
+			Log(Log::DEBUG)<<__func__<<": "<<fileNames[i]<<" as YAML image";
+			tmpImage = openImageYml(fileNames[i]);
+		}
+		else
+		{
+			Log(Log::DEBUG)<<__func__<<": "<<fileNames[i]<<" as png image";
+			tmpImage = openImageImg(fileNames[i]);
+		}
 		if(tmpImage.data)
 			images.push_back(tmpImage);
 		else 
-			Log(Log::WARN)<<"can not read image "<<i<<" from "<<fileNames[i]<<'\n';
+			Log(Log::WARN)<<"can not read image "<<i<<" from "<<fileNames[i];
 	}
 	return images;
 }
@@ -129,7 +138,7 @@ std::vector<cv::Mat> loadImages(char** fileNames)
 int perfromOperation(int operation, char** fileNames, const Config& config)
 {
 	std::vector<cv::Mat> inImages;
-	if(operation == CREATE_MAP || operation == APPLY_MAP || operation == APPLY_CURVE)
+	if(operation == CREATE_MAP || operation == APPLY_MAP || operation == APPLY_CURVE || operation == SHOW_IMAGE)
 	{
 		inImages = loadImages(fileNames);
 		
@@ -323,6 +332,16 @@ int perfromOperation(int operation, char** fileNames, const Config& config)
 		
 		std::cout<<"Curve saved to "<<(!config.output.empty() ? config.output : "curve.mat")<<'\n';
 	}
+	else if(operation == SHOW_IMAGE)
+	{
+		cv::namedWindow("Show Image", cv::WINDOW_NORMAL );
+		for(size_t i = 0; i < inImages.size(); ++i)
+		{
+			cv::imshow("Show Image", inImages[i]);
+			cv::waitKey(0);
+		}
+		cv::destroyWindow("Show Image"); 
+	}
 	else if(operation == EXIT)
 		return -1;
 	return 0;
diff --git a/src/matutils.cpp b/src/matutils.cpp
index bf4d893..5bea21d 100644
--- a/src/matutils.cpp
+++ b/src/matutils.cpp
@@ -178,7 +178,7 @@ bool findClosest(size_t& index, const cv::Point2f point, const std::vector<cv::P
 	return found;
 }
 
-void interpolateMissing(cv::Mat& mat)
+void interpolateMissingOnX(cv::Mat& mat)
 {
 	assert(mat.type() == CV_32FC1);
 	
@@ -223,35 +223,48 @@ void interpolateMissing(cv::Mat& mat)
 	}
 }
 
-void fillMissing(cv::Mat& mat)
+void interpolateMissingOnY(cv::Mat& mat)
 {
 	assert(mat.type() == CV_32FC1);
 	
-	bool finished = true;
-	for(int y = 0; y < mat.rows; y++)
+	for(int x = 0; x < mat.cols; ++x)
 	{
-		float* col = mat.ptr<float>(y);
-		for(int x = 0; x < mat.cols; ++x)
+		for(int y = 0; y < mat.rows; y++)
 		{
-			if(col[x] < 0 && col[x] > -2) 
+			if(mat.at<float>(y,x) < 0) 
 			{
-				if(y > 0 && mat.at<float>(y-1,x) >= 0)
+				int closestA = -1;
+				int closestB = -1;
+				int dist = std::numeric_limits<int>::max();
+				for(int i = 0; i < mat.rows; i++)
 				{
-					col[x] = mat.at<float>(y-1,x);
-					finished = false;
+					if(i != closestA && mat.at<float>(i,x) >= 0 && abs(i-y) <= dist)
+					{
+						closestB = closestA;
+						closestA = i;
+						dist = abs(i-y);
+					}
 				}
-				else if(y < mat.rows-1 && mat.at<float>(y+1,x) >= 0)
+				if(closestA < 0 || closestB < 0)
 				{
-					col[x] = mat.at<float>(y+1,x);
-					finished = false;
+					closestA = -1;
+					closestB = -1;
+					dist = std::numeric_limits<int>::max();
+					for(int i = mat.rows-1; i >= 0; --i)
+					{
+						if(i != closestA && mat.at<float>(i,x) >= 0 && abs(i-y) <= dist)
+						{
+							closestB = closestA;
+							closestA = i;
+							dist = abs(i-y);
+						}
+					}
 				}
-				if(col[x] > 0 && ((x+1 < mat.cols && col[x] > col[x+1]) || (x > 0 && col[x] < col[x-1])))
-					col[x] = -2;
+				float slope = (mat.at<float>(closestB,x) - mat.at<float>(closestA,x))/(closestB-closestA);
+				mat.at<float>(y,x) = mat.at<float>(closestA,x) - (closestA-y)*slope;
 			}
 		}
 	}
-	if(!finished) 
-		fillMissing(mat);
 }
 
 bool findDeadSpace(const cv::Mat& mat, cv::Rect& roi)
@@ -295,8 +308,8 @@ void removeSparseCols(cv::Mat& mat, float reject, bool front)
 	int count = 0;
 	for(int y = 0; y < mat.rows; ++y)
 	{
-		if(front && mat.at<float>(y,0) >= 0) ++count;
-		else if(mat.at<float>(y,mat.cols-1) >= 0) ++count;
+		if((front && mat.at<float>(y,0) >= 0) || (!front && mat.at<float>(y,mat.cols-1) >= 0))
+			++count;
 	}
 	cv::Rect roi;
 	bool rej = (count < mat.rows*reject);
diff --git a/src/matutils.h b/src/matutils.h
index 3d711e9..6020967 100644
--- a/src/matutils.h
+++ b/src/matutils.h
@@ -37,9 +37,9 @@ bool findClosest(size_t& index, const cv::Point2f point, const std::vector<cv::P
 
 void thompsonTauTest(const std::vector<float>& in, std::vector<size_t>& outliers, float criticalValue);
 
-void interpolateMissing(cv::Mat& mat);
+void interpolateMissingOnX(cv::Mat& mat);
 
-void fillMissing(cv::Mat& mat);
+void interpolateMissingOnY(cv::Mat& mat);
 
 bool findDeadSpace(const cv::Mat& mat, cv::Rect& roi);
 
diff --git a/src/unwrap.cpp b/src/unwrap.cpp
index 7414c96..d0bfac4 100644
--- a/src/unwrap.cpp
+++ b/src/unwrap.cpp
@@ -66,7 +66,7 @@ static void sortIntoRemapMaps(const std::vector<DetectedPoint>& points, cv::Mat&
 	xMat = -1;
 	yMat = -1;
 	
-	Log(Log::DEBUG)<<"Grid: "<<xGridSize<<'x'<<yGridSize;
+	Log(Log::DEBUG)<<__func__<<"Grid: "<<xGridSize<<'x'<<yGridSize;
 	
 	for(int y = 0; y < xMat.rows; y++)
 	{
@@ -103,16 +103,14 @@ bool createRemapMap(const cv::Mat& image, RemapMap& out, const std::vector<Detec
 	removeSparseCols(out.yMat, 0.75, true);
 	removeSparseCols(out.yMat, 0.75, false);
 	Log(Log::DEBUG)<<__func__<<": xMat rejcted\n"<<out.xMat;
-	fillMissing(out.xMat);
 	Log(Log::DEBUG)<<__func__<<": xMat filled\n"<<out.xMat;
-	interpolateMissing(out.xMat);
-	interpolateMissing(out.yMat);
+	interpolateMissingOnY(out.xMat);
+	interpolateMissingOnY(out.yMat);
 
 	sanityCheckMap(out.xMat, 0, image.cols-1, -1, -1);
 	sanityCheckMap(out.yMat, 0, image.rows-1, -1, -1);
-	fillMissing(out.xMat);
-	interpolateMissing(out.xMat);
-	interpolateMissing(out.yMat);
+	interpolateMissingOnY(out.xMat);
+	interpolateMissingOnY(out.yMat);
 	sanityCheckMap(out.xMat, 0, image.cols-1, 0, image.cols-1);
 	sanityCheckMap(out.yMat, 0, image.rows-1, 0, image.rows-1);
 	
@@ -198,10 +196,10 @@ cv::Mat simpleStich(const std::vector<RemapedImage>& images)
 		if(topLeft.y > image.origin.y)
 			topLeft.y = image.origin.y;
 		
-		Log(Log::DEBUG)<<"image: "<<image.image.rows<<'x'<<image.image.cols<<" at "<<image.origin.x<<'x'<<image.origin.y;
+		Log(Log::DEBUG)<<__func__<<"image: "<<image.image.rows<<'x'<<image.image.cols<<" at "<<image.origin.x<<'x'<<image.origin.y;
 	}
 	
-	Log(Log::DEBUG)<<"outputSize: "<<outputSize;
+	Log(Log::DEBUG)<<__func__<<"outputSize: "<<outputSize;
 	
 	cv::Mat out(outputSize, images[0].image.type(), cv::Scalar::all(0));
 	
@@ -219,9 +217,13 @@ cv::Mat simpleStich(const std::vector<RemapedImage>& images)
 
 cv::Mat stich(std::vector<RemapedImage>& images, bool seamAdjust)
 {
+	assert(images.size() > 0);
+
 	for(auto& image : images)
 		assert(image.image.type() == CV_8UC3 || image.image.type() == CV_8UC1);
 
+	Log(Log::DEBUG)<<__func__<<" got "<<images.size()<<" of type "<<images[0].image.type()<<" size "<<images[0].image.size();
+
 	for(auto& image : images)
 	{
 		if(image.image.type() == CV_8UC1)