 		)
 		{
 			
-			// TODO: Maybe add special cases if only one or two neighbors are missing?			
-
-			switch (sourceInfo.interpolation)
+			//targetInfo.interpolation = himan::kNearestPointValue;
+
+			switch (targetInfo.interpolation)
 			{
 				case himan::kBiLinear:
-				{
-					// Find all four neighboring points
-
-					point a(floor(gp.x), ceil(gp.y));
-					point b(ceil(gp.x), ceil(gp.y));
-					point c(floor(gp.x), floor(gp.y));
-					point d(ceil(gp.x), floor(gp.y));
-
-					// Neighbor values
-
-					double av = d_source[Index(a,sourceInfo.size_x)];
-					double bv = d_source[Index(b,sourceInfo.size_x)];
-					double cv = d_source[Index(c,sourceInfo.size_x)];
-					double dv = d_source[Index(d,sourceInfo.size_x)];	
-
-					// Distance of interpolated point to neighboring points
-
-					point dist(gp.x - c.x, gp.y - c.y);
-
-					if (av != kFloatMissing && bv != kFloatMissing && cv != kFloatMissing && dv != kFloatMissing)
-					{
-						interp = BiLinear(dist.x, dist.y, av, bv, cv, dv);
-					}
-
-#ifdef EXTRADEBUG
-					// Neighbor point indexes in linear format
-
-					int aidx = Index(a,sourceInfo.size_x);
-					int bidx = Index(b,sourceInfo.size_x);
-					int cidx = Index(c,sourceInfo.size_x);
-					int didx = Index(d,sourceInfo.size_x);
-
-					if (i == 0 && j == 0)
-					{
-						printf("x:%d y:%d gpx:%f gpy:%f\n", i, j, gp.x, gp.y);
-						printf("a x:%d y:%d val:%f\n", int(a.x), int(a.y), av);
-						printf("b x:%d y:%d val:%f\n", int(b.x), int(b.y), bv);
-						printf("c x:%d y:%d val:%f\n", int(c.x), int(c.y), cv);
-						printf("d x:%d y:%d val:%f\n", int(d.x), int(d.y), dv);
-						printf("dist x:%f y:%f\n", dist.x, dist.y);
-						printf("interp:%f\n", interp);
-					}
-#endif
-				}
+					interp = BiLinearInterpolation(d_source, sourceInfo, gp);
 					break;
 					
 				case himan::kNearestPoint:
-					int rx = rint(gp.x);
-					int ry = rint(gp.y);
-					
-					assert(rx >= 0 && rx < sourceInfo.size_x);
-					assert(ry >= 0 && ry < sourceInfo.size_y);
-					
-					interp = d_source[Index(rx,ry,sourceInfo.size_x)];
-					
+					interp = NearestPointInterpolation(d_source, sourceInfo, gp);
 					break;
 					
+				case himan::kNearestPointValue:
+					interp = NearestPointValueInterpolation(d_source, sourceInfo, gp);
+					break;
 			}
 		}
 

}

__device__
double Linear(double dx, double left, double right)
{
	// return (1 - dx) * left + dx * right;

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

const double kEpsilon = 1e-5;

struct point
{
