 
 	double ret = kFloatMissing;
 
+	// All values present, regular bilinear interpolation
+	
 	if (av != kFloatMissing && bv != kFloatMissing && cv != kFloatMissing && dv != kFloatMissing)
 	{
 		ret = BiLinear(dist.x, dist.y, av, bv, cv, dv);
 	}
 
-	// These "triangulation" methods have been copied from NFmiInterpolation.cpp
+	// x or y is at grid edge
+	
+	else if (fabs(dist.y) < kEpsilon && cv != kFloatMissing && dv != kFloatMissing)
+	{
+		ret = Linear(dist.x, cv, dv);
+	}
+
+	else if (fabs(dist.y - 1) < kEpsilon && av != kFloatMissing && bv != kFloatMissing)
+	{
+		ret = Linear(dist.x, av, bv);
+	}
+	
+	else if (fabs(dist.x) < kEpsilon && cv != kFloatMissing && av != kFloatMissing)
+	{
+		ret = Linear(dist.y, cv, dv);
+	}
+
+	else if (fabs(dist.x - 1) < kEpsilon && av != kFloatMissing && bv != kFloatMissing)
+	{
+		ret = Linear(dist.y, dv, bv);
+	}
+		
+	// One point missing; these "triangulation" methods have been copied from NFmiInterpolation.cpp
 
 	else if (av == kFloatMissing && bv != kFloatMissing && cv != kFloatMissing && dv != kFloatMissing)
 	{

			// if interpolated grid points are negative, it means that we are outside
			// of the source area
				
			gp.x >= 0 && gp.y >= 0 &&
				
			// if interpolated grid points are larger than source grid in x or y
			// direction, it means again that we are outside of the area
			
			gp.x < (sourceInfo.size_x-1) && gp.y < (sourceInfo.size_y-1)
		)
		{
			

#include "cuda_helper.h"
#include <thrust/sort.h>

const double kEpsilon = 1e-6;

struct point
{
