X-Git-Url: http://git.xonotic.org/?p=xonotic%2Fxonotic.git;a=blobdiff_plain;f=misc%2Ftools%2Ffft-normalmap-to-heightmap.c;h=88ac27c9fde1e4288d3907f198669c0b06a6491a;hp=c133515300de23ef9ec95cd18bba010f6ec11a57;hb=395a8b52d5e0b3b6a102cb7e2383902a8c3f0ae2;hpb=c151bbe67628c30839d1ba0ce518a610e626770a

diff --git a/misc/tools/fft-normalmap-to-heightmap.c b/misc/tools/fft-normalmap-to-heightmap.c
index c1335153..88ac27c9 100644
--- a/misc/tools/fft-normalmap-to-heightmap.c
+++ b/misc/tools/fft-normalmap-to-heightmap.c
@@ -33,9 +33,14 @@
 
 #include <fftw3.h>
 
-void nmap_to_hmap(unsigned char *map, int w, int h, double scale, double offset)
+#define TWO_PI (4*atan2(1,1) * 2)
+
+void nmap_to_hmap(unsigned char *map, const unsigned char *refmap, int w, int h, double scale, double offset, const double *filter, int filterw, int filterh)
 {
 	int x, y;
+	int i, j;
+	double fx, fy;
+	double ffx, ffy;
 	double nx, ny, nz;
 	double v, vmin, vmax;
 #ifndef C99
@@ -46,9 +51,9 @@ void nmap_to_hmap(unsigned char *map, int w, int h, double scale, double offset)
 	fftw_complex *imgspace2 = fftw_malloc(w*h * sizeof(fftw_complex));
 	fftw_complex *freqspace1 = fftw_malloc(w*h * sizeof(fftw_complex));
 	fftw_complex *freqspace2 = fftw_malloc(w*h * sizeof(fftw_complex));
-	fftw_plan i12f1 = fftw_plan_dft_2d(w, h, imgspace1, freqspace1, FFTW_FORWARD, FFTW_ESTIMATE);
-	fftw_plan i22f2 = fftw_plan_dft_2d(w, h, imgspace2, freqspace2, FFTW_FORWARD, FFTW_ESTIMATE);
-	fftw_plan f12i1 = fftw_plan_dft_2d(w, h, freqspace1, imgspace1, FFTW_BACKWARD, FFTW_ESTIMATE);
+	fftw_plan i12f1 = fftw_plan_dft_2d(h, w, imgspace1, freqspace1, FFTW_FORWARD, FFTW_ESTIMATE);
+	fftw_plan i22f2 = fftw_plan_dft_2d(h, w, imgspace2, freqspace2, FFTW_FORWARD, FFTW_ESTIMATE);
+	fftw_plan f12i1 = fftw_plan_dft_2d(h, w, freqspace1, imgspace1, FFTW_BACKWARD, FFTW_ESTIMATE);
 
 	for(y = 0; y < h; ++y)
 	for(x = 0; x < w; ++x)
@@ -60,18 +65,18 @@ void nmap_to_hmap(unsigned char *map, int w, int h, double scale, double offset)
 		 * n_z = -dh/dh = -1
 		 * BUT: darkplaces uses inverted normals, n_y actually is dh/dy by image pixel coordinates
 		 */
-		nx = (int)map[(w*y+x)*4+2] - 127.5;
-		ny = (int)map[(w*y+x)*4+1] - 127.5;
-		nz = (int)map[(w*y+x)*4+0] - 127.5;
+		nx = ((int)map[(w*y+x)*4+2] - 127.5) / 128;
+		ny = ((int)map[(w*y+x)*4+1] - 127.5) / 128;
+		nz = ((int)map[(w*y+x)*4+0] - 127.5) / 128;
 
 		/* reconstruct the derivatives from here */
 #ifdef C99
-		imgspace1[(w*y+x)] =  nx / nz; /* = dz/dx */
-		imgspace2[(w*y+x)] = -ny / nz; /* = dz/dy */
+		imgspace1[(w*y+x)] =  nx / nz * w; /* = dz/dx */
+		imgspace2[(w*y+x)] = -ny / nz * h; /* = dz/dy */
 #else
-		imgspace1[(w*y+x)][0] =  nx / nz; /* = dz/dx */
+		imgspace1[(w*y+x)][0] =  nx / nz * w; /* = dz/dx */
 		imgspace1[(w*y+x)][1] = 0;
-		imgspace2[(w*y+x)][0] = -ny / nz; /* = dz/dy */
+		imgspace2[(w*y+x)][0] = -ny / nz * h; /* = dz/dy */
 		imgspace2[(w*y+x)][1] = 0;
 #endif
 	}
@@ -84,69 +89,183 @@ void nmap_to_hmap(unsigned char *map, int w, int h, double scale, double offset)
 	for(y = 0; y < h; ++y)
 	for(x = 0; x < w; ++x)
 	{
-		int fx = x;
-		int fy = y;
-		if(fx > w/2)
-			fx -= w;
-		if(fy > h/2)
-			fy -= h;
+		fx = x * 1.0 / w;
+		fy = y * 1.0 / h;
+		if(filter)
+		{
+			// discontinous case
+			// we must invert whatever "filter" would do on (x, y)!
 #ifdef C99
-		if(fx||fy)
-			freqspace1[(w*y+x)] = I * (fx * freqspace1[(w*y+x)] + fy * freqspace2[(w*y+x)]) / (fx*fx + fy*fy);
-		else
-			freqspace1[(w*y+x)] = 0;
+			fftw_complex response_x = 0;
+			fftw_complex response_y = 0;
+			double sum;
+			for(i = -filterh / 2; i <= filterh / 2; ++i)
+				for(j = -filterw / 2; j <= filterw / 2; ++j)
+				{
+					response_x += filter[(i + filterh / 2) * filterw + j + filterw / 2] * cexp(-_Complex_I * TWO_PI * (j * fx + i * fy));
+					response_y += filter[(i + filterh / 2) * filterw + j + filterw / 2] * cexp(-_Complex_I * TWO_PI * (i * fx + j * fy));
+				}
+
+			// we know:
+			//   fourier(df/dx)_xy = fourier(f)_xy * response_x
+			//   fourier(df/dy)_xy = fourier(f)_xy * response_y
+			// mult by conjugate of response_x, response_y:
+			//   conj(response_x) * fourier(df/dx)_xy = fourier(f)_xy * |response_x^2|
+			//   conj(response_y) * fourier(df/dy)_xy = fourier(f)_xy * |response_y^2|
+			// and
+			//   fourier(f)_xy = (conj(response_x) * fourier(df/dx)_xy + conj(response_y) * fourier(df/dy)_xy) / (|response_x|^2 + |response_y|^2)
+
+			sum = cabs(response_x) * cabs(response_x) + cabs(response_y) * cabs(response_y);
+
+			if(sum > 0)
+				freqspace1[(w*y+x)] = (conj(response_x) * freqspace1[(w*y+x)] + conj(response_y) * freqspace2[(w*y+x)]) / sum;
+			else
+				freqspace1[(w*y+x)] = 0;
 #else
-		if(fx||fy)
-		{
-			save = freqspace1[(w*y+x)][0];
-			freqspace1[(w*y+x)][0] = -(fx * freqspace1[(w*y+x)][1] + fy * freqspace2[(w*y+x)][1]) / (fx*fx + fy*fy);
-			freqspace1[(w*y+x)][1] = (fx * save + fy * freqspace2[(w*y+x)][0]) / (fx*fx + fy*fy);
+			fftw_complex response_x = {0, 0};
+			fftw_complex response_y = {0, 0};
+			double sum;
+			for(i = -filterh / 2; i <= filterh / 2; ++i)
+				for(j = -filterw / 2; j <= filterw / 2; ++j)
+				{
+					response_x[0] += filter[(i + filterh / 2) * filterw + j + filterw / 2] * cos(-TWO_PI * (j * fx + i * fy));
+					response_x[1] += filter[(i + filterh / 2) * filterw + j + filterw / 2] * sin(-TWO_PI * (j * fx + i * fy));
+					response_y[0] += filter[(i + filterh / 2) * filterw + j + filterw / 2] * cos(-TWO_PI * (i * fx + j * fy));
+					response_y[1] += filter[(i + filterh / 2) * filterw + j + filterw / 2] * sin(-TWO_PI * (i * fx + j * fy));
+				}
+
+			sum = response_x[0] * response_x[0] + response_x[1] * response_x[1]
+			    + response_y[0] * response_y[0] + response_y[1] * response_y[1];
+
+			if(sum > 0)
+			{
+				double s = freqspace1[(w*y+x)][0];
+				freqspace1[(w*y+x)][0] = (response_x[0] * s                      + response_x[1] * freqspace1[(w*y+x)][1] + response_y[0] * freqspace2[(w*y+x)][0] + response_y[1] * freqspace2[(w*y+x)][1]) / sum;
+				freqspace1[(w*y+x)][1] = (response_x[0] * freqspace1[(w*y+x)][1] - response_x[1] * s                      + response_y[0] * freqspace2[(w*y+x)][1] - response_y[1] * freqspace2[(w*y+x)][0]) / sum;
+			}
+			else
+			{
+				freqspace1[(w*y+x)][0] = 0;
+				freqspace1[(w*y+x)][1] = 0;
+			}
+#endif
 		}
 		else
 		{
-			freqspace1[(w*y+x)][0] = 0;
-			freqspace1[(w*y+x)][1] = 0;
-		}
+			// continuous integration case
+			if(fx > 0.5)
+				fx -= 1;
+			if(fy > 0.5)
+				fy -= 1;
+			/* these must have the same sign as fx and fy (so ffx*fx + ffy*fy is nonzero), otherwise do not matter */
+			/* it basically decides how artifacts are distributed */
+			ffx = fx;
+			ffy = fy;
+#ifdef C99
+			if(fx||fy)
+				freqspace1[(w*y+x)] = _Complex_I * (ffx * freqspace1[(w*y+x)] + ffy * freqspace2[(w*y+x)]) / (ffx*fx + ffy*fy) / TWO_PI;
+			else
+				freqspace1[(w*y+x)] = 0;
+#else
+			if(fx||fy)
+			{
+				save = freqspace1[(w*y+x)][0];
+				freqspace1[(w*y+x)][0] = -(ffx * freqspace1[(w*y+x)][1] + ffy * freqspace2[(w*y+x)][1]) / (ffx*fx + ffy*fy) / TWO_PI;
+				freqspace1[(w*y+x)][1] =  (ffx * save + ffy * freqspace2[(w*y+x)][0]) / (ffx*fx + ffy*fy) / TWO_PI;
+			}
+			else
+			{
+				freqspace1[(w*y+x)][0] = 0;
+				freqspace1[(w*y+x)][1] = 0;
+			}
 #endif
+		}
 	}
 
 	fftw_execute(f12i1);
 
-	if(scale == 0)
+	/* renormalize, find min/max */
+	vmin = vmax = 0;
+	for(y = 0; y < h; ++y)
+	for(x = 0; x < w; ++x)
 	{
 #ifdef C99
-		vmin = vmax = creal(imgspace1[0]);
+		v = creal(imgspace1[(w*y+x)] /= pow(w*h, 1.5));
 #else
-		vmin = vmax = imgspace1[0][0];
+		v = (imgspace1[(w*y+x)][0] /= pow(w*h, 1.5));
+		// imgspace1[(w*y+x)][1] /= pow(w*h, 1.5);
+		// this value is never used
 #endif
+		if(v < vmin || (x == 0 && y == 0))
+			vmin = v;
+		if(v > vmax || (x == 0 && y == 0))
+			vmax = v;
+	}
+
+	if(refmap)
+	{
+		double f, a;
+		double o, s;
+		double sa, sfa, sffa, sfva, sva;
+		double mi, ma;
+		sa = sfa = sffa = sfva = sva = 0;
+		mi = 1;
+		ma = -1;
 		for(y = 0; y < h; ++y)
 		for(x = 0; x < w; ++x)
 		{
+			a = (int)refmap[(w*y+x)*4+3];
+			v = (refmap[(w*y+x)*4+0]*0.114 + refmap[(w*y+x)*4+1]*0.587 + refmap[(w*y+x)*4+2]*0.299);
+			v = (v - 128.0) / 127.0;
 #ifdef C99
-			v = creal(imgspace1[(w*y+x)]);
+			f = creal(imgspace1[(w*y+x)]);
 #else
-			v = imgspace1[(w*y+x)][0];
+			f = imgspace1[(w*y+x)][0];
 #endif
-			if(v < vmin)
-				vmin = v;
-			if(v > vmax)
-				vmax = v;
+			if(a <= 0)
+				continue;
+			if(v < mi)
+				mi = v;
+			if(v > ma)
+				ma = v;
+			sa += a;
+			sfa += f*a;
+			sffa += f*f*a;
+			sfva += f*v*a;
+			sva += v*a;
+		}
+		if(mi < ma)
+		{
+			/* linear regression ftw */
+			o = (sfa*sfva - sffa*sva) / (sfa*sfa-sa*sffa);
+			s = (sfa*sva - sa*sfva) / (sfa*sfa-sa*sffa);
+		}
+		else /* all values of v are equal, so we cannot get scale; we can still get offset */
+		{
+			o = ((sva - sfa) / sa);
+			s = 1;
 		}
-
-		vmin /= (w*h);
-		vmax /= (w*h);
 
+		/*
+		 * now apply user-given offset and scale to these values
+		 * (x * s + o) * scale + offset
+		 * x * s * scale + o * scale + offset
+		 */
+		offset += o * scale;
+		scale *= s;
+	}
+	else if(scale == 0)
+	{
 		/*
 		 * map vmin to -1
 		 * map vmax to +1
 		 */
 		scale = 2 / (vmax - vmin);
 		offset = -(vmax + vmin) / (vmax - vmin);
-
-		printf("Autocomputed scale: %f\nAutocomputed offset: %f\n", scale, offset);
 	}
 
-	scale /= (w*h);
+	printf("Min: %f\nAvg: %f\nMax: %f\nScale: %f\nOffset: %f\nScaled-Min: %f\nScaled-Avg: %f\nScaled-Max: %f\n", 
+		vmin, 0.0, vmax, scale, offset, vmin * scale + offset, offset, vmax * scale + offset);
 
 	for(y = 0; y < h; ++y)
 	for(x = 0; x < w; ++x)
@@ -161,12 +280,6 @@ void nmap_to_hmap(unsigned char *map, int w, int h, double scale, double offset)
 			v = -1;
 		if(v > 1)
 			v = 1;
-		/*
-		map[(w*y+x)*4+0] = 128 + 127 * v;
-		map[(w*y+x)*4+1] = 128 + 127 * v;
-		map[(w*y+x)*4+2] = 128 + 127 * v;
-		map[(w*y+x)*4+3] = 255;
-		*/
 		map[(w*y+x)*4+3] = floor(128.5 + 127 * v);
 	}
 
@@ -180,6 +293,192 @@ void nmap_to_hmap(unsigned char *map, int w, int h, double scale, double offset)
 	fftw_free(imgspace1);
 }
 
+void hmap_to_nmap(unsigned char *map, int w, int h, int src_chan, double scale)
+{
+	int x, y;
+	double fx, fy;
+	double nx, ny, nz;
+	double v;
+#ifndef C99
+	double save;
+#endif
+
+	fftw_complex *imgspace1 = fftw_malloc(w*h * sizeof(fftw_complex));
+	fftw_complex *imgspace2 = fftw_malloc(w*h * sizeof(fftw_complex));
+	fftw_complex *freqspace1 = fftw_malloc(w*h * sizeof(fftw_complex));
+	fftw_complex *freqspace2 = fftw_malloc(w*h * sizeof(fftw_complex));
+	fftw_plan i12f1 = fftw_plan_dft_2d(h, w, imgspace1, freqspace1, FFTW_FORWARD, FFTW_ESTIMATE);
+	fftw_plan f12i1 = fftw_plan_dft_2d(h, w, freqspace1, imgspace1, FFTW_BACKWARD, FFTW_ESTIMATE);
+	fftw_plan f22i2 = fftw_plan_dft_2d(h, w, freqspace2, imgspace2, FFTW_BACKWARD, FFTW_ESTIMATE);
+
+	for(y = 0; y < h; ++y)
+	for(x = 0; x < w; ++x)
+	{
+		switch(src_chan)
+		{
+			case 0:
+			case 1:
+			case 2:
+			case 3:
+				v = map[(w*y+x)*4+src_chan];
+				break;
+			case 4:
+				v = (map[(w*y+x)*4+0] + map[(w*y+x)*4+1] + map[(w*y+x)*4+2]) / 3;
+				break;
+			default:
+			case 5:
+				v = (map[(w*y+x)*4+0]*0.114 + map[(w*y+x)*4+1]*0.587 + map[(w*y+x)*4+2]*0.299);
+				break;
+		}
+#ifdef C99
+		imgspace1[(w*y+x)] = (v - 128.0) / 127.0;
+#else
+		imgspace1[(w*y+x)][0] = (v - 128.0) / 127.0;
+		imgspace1[(w*y+x)][1] = 0;
+#endif
+		if(v < 1)
+			v = 1; /* do not write alpha zero */
+		map[(w*y+x)*4+3] = floor(v + 0.5);
+	}
+
+	/* see http://www.gamedev.net/community/forums/topic.asp?topic_id=561430 */
+
+	fftw_execute(i12f1);
+	
+	for(y = 0; y < h; ++y)
+	for(x = 0; x < w; ++x)
+	{
+		fx = x;
+		fy = y;
+		if(fx > w/2)
+			fx -= w;
+		if(fy > h/2)
+			fy -= h;
+#ifdef DISCONTINUOUS
+		fx = sin(fx * TWO_PI / w);
+		fy = sin(fy * TWO_PI / h);
+#else
+#ifdef C99
+		/* a lowpass to prevent the worst */
+		freqspace1[(w*y+x)] *= 1 - pow(abs(fx) / (double)(w/2), 1);
+		freqspace1[(w*y+x)] *= 1 - pow(abs(fy) / (double)(h/2), 1);
+#else
+		/* a lowpass to prevent the worst */
+		freqspace1[(w*y+x)][0] *= 1 - pow(abs(fx) / (double)(w/2), 1);
+		freqspace1[(w*y+x)][1] *= 1 - pow(abs(fx) / (double)(w/2), 1);
+		freqspace1[(w*y+x)][0] *= 1 - pow(abs(fy) / (double)(h/2), 1);
+		freqspace1[(w*y+x)][1] *= 1 - pow(abs(fy) / (double)(h/2), 1);
+#endif
+#endif
+#ifdef C99
+		freqspace2[(w*y+x)] = TWO_PI*_Complex_I * fy * freqspace1[(w*y+x)]; /* y derivative */
+		freqspace1[(w*y+x)] = TWO_PI*_Complex_I * fx * freqspace1[(w*y+x)]; /* x derivative */
+#else
+		freqspace2[(w*y+x)][0] = -TWO_PI * fy * freqspace1[(w*y+x)][1]; /* y derivative */
+		freqspace2[(w*y+x)][1] =  TWO_PI * fy * freqspace1[(w*y+x)][0];
+		save = freqspace1[(w*y+x)][0];
+		freqspace1[(w*y+x)][0] = -TWO_PI * fx * freqspace1[(w*y+x)][1]; /* x derivative */
+		freqspace1[(w*y+x)][1] =  TWO_PI * fx * save;
+#endif
+	}
+
+	fftw_execute(f12i1);
+	fftw_execute(f22i2);
+
+	scale /= (w*h);
+
+	for(y = 0; y < h; ++y)
+	for(x = 0; x < w; ++x)
+	{
+#ifdef C99
+		nx = creal(imgspace1[(w*y+x)]);
+		ny = creal(imgspace2[(w*y+x)]);
+#else
+		nx = imgspace1[(w*y+x)][0];
+		ny = imgspace2[(w*y+x)][0];
+#endif
+		nx /= w;
+		ny /= h;
+		nz = -1 / scale;
+		v = -sqrt(nx*nx + ny*ny + nz*nz);
+		nx /= v;
+		ny /= v;
+		nz /= v;
+		ny = -ny; /* DP inverted normals */
+		map[(w*y+x)*4+2] = floor(128 + 127.5 * nx);
+		map[(w*y+x)*4+1] = floor(128 + 127.5 * ny);
+		map[(w*y+x)*4+0] = floor(128 + 127.5 * nz);
+	}
+
+	fftw_destroy_plan(i12f1);
+	fftw_destroy_plan(f12i1);
+	fftw_destroy_plan(f22i2);
+
+	fftw_free(freqspace2);
+	fftw_free(freqspace1);
+	fftw_free(imgspace2);
+	fftw_free(imgspace1);
+}
+
+void hmap_to_nmap_local(unsigned char *map, int w, int h, int src_chan, double scale, const double *filter, int filterw, int filterh)
+{
+	int x, y;
+	double nx, ny, nz;
+	double v;
+	int i, j;
+	double *img_reduced = malloc(w*h * sizeof(double));
+
+	for(y = 0; y < h; ++y)
+	for(x = 0; x < w; ++x)
+	{
+		switch(src_chan)
+		{
+			case 0:
+			case 1:
+			case 2:
+			case 3:
+				v = map[(w*y+x)*4+src_chan];
+				break;
+			case 4:
+				v = (map[(w*y+x)*4+0] + map[(w*y+x)*4+1] + map[(w*y+x)*4+2]) / 3;
+				break;
+			default:
+			case 5:
+				v = (map[(w*y+x)*4+0]*0.114 + map[(w*y+x)*4+1]*0.587 + map[(w*y+x)*4+2]*0.299);
+				break;
+		}
+		img_reduced[(w*y+x)] = (v - 128.0) / 127.0;
+		if(v < 1)
+			v = 1; /* do not write alpha zero */
+		map[(w*y+x)*4+3] = floor(v + 0.5);
+	}
+
+	for(y = 0; y < h; ++y)
+	for(x = 0; x < w; ++x)
+	{
+		nz = -1 / scale;
+		nx = ny = 0;
+
+		for(i = -filterh / 2; i <= filterh / 2; ++i)
+			for(j = -filterw / 2; j <= filterw / 2; ++j)
+			{
+				nx += img_reduced[w*((y+i+h)%h)+(x+j+w)%w] * filter[(i + filterh / 2) * filterw + j + filterw / 2];
+				ny += img_reduced[w*((y+j+h)%h)+(x+i+w)%w] * filter[(i + filterh / 2) * filterw + j + filterw / 2];
+			}
+
+		v = -sqrt(nx*nx + ny*ny + nz*nz);
+		nx /= v;
+		ny /= v;
+		nz /= v;
+		ny = -ny; /* DP inverted normals */
+		map[(w*y+x)*4+2] = floor(128 + 127.5 * nx);
+		map[(w*y+x)*4+1] = floor(128 + 127.5 * ny);
+		map[(w*y+x)*4+0] = floor(128 + 127.5 * nz);
+	}
+
+	free(img_reduced);
+}
+
 unsigned char *FS_LoadFile(const char *fn, int *len)
 {
 	unsigned char *buf = NULL;
@@ -663,16 +962,74 @@ int Image_WriteTGABGRA (const char *filename, int width, int height, const unsig
 
 int usage(const char *me)
 {
-	printf("Usage: %s <infile.tga> <outfile.tga> [<scale> [<offset>]]\n", me);
+	printf("Usage: %s <infile_norm.tga> <outfile_normandheight.tga> filtertype [<scale> [<offset> [<infile_ref.tga>]]] (get heightmap from normalmap)\n", me);
+	printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> filtertype -1 [<scale>] (read from B)\n", me);
+	printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> filtertype -2 [<scale>] (read from G)\n", me);
+	printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> filtertype -3 [<scale>] (read from R)\n", me);
+	printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> filtertype -4 [<scale>] (read from A)\n", me);
+	printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> filtertype -5 [<scale>] (read from (R+G+B)/3)\n", me);
+	printf("or:    %s <infile_height.tga> <outfile_normandheight.tga> filtertype -6 [<scale>] (read from Y)\n", me);
 	return 1;
 }
 
+static const double filter_scharr3[3][3] = {
+	{  -3/32.0, 0,  3/32.0 },
+	{ -10/32.0, 0, 10/32.0 },
+	{  -3/32.0, 0,  3/32.0 }
+};
+
+static const double filter_prewitt3[3][3] = {
+	{ -1/6.0, 0, 1/6.0 },
+	{ -1/6.0, 0, 1/6.0 },
+	{ -1/6.0, 0, 1/6.0 }
+};
+
+// pathologic for inverting
+static const double filter_sobel3[3][3] = {
+	{ -1/8.0, 0, 1/8.0 },
+	{ -2/8.0, 0, 2/8.0 },
+	{ -1/8.0, 0, 1/8.0 }
+};
+
+// pathologic for inverting
+static const double filter_sobel5[5][5] = {
+	{ -1/128.0,  -2/128.0, 0,  2/128.0, 1/128.0 },
+	{ -4/128.0,  -8/128.0, 0,  8/128.0, 4/128.0 },
+	{ -6/128.0, -12/128.0, 0, 12/128.0, 6/128.0 },
+	{ -4/128.0,  -8/128.0, 0,  8/128.0, 4/128.0 },
+	{ -1/128.0,  -2/128.0, 0,  2/128.0, 1/128.0 }
+};
+
+// pathologic for inverting
+static const double filter_prewitt5[5][5] = {
+	{ -1/40.0, -2/40.0, 0, 2/40.0, 1/40.0 },
+	{ -1/40.0, -2/40.0, 0, 2/40.0, 1/40.0 },
+	{ -1/40.0, -2/40.0, 0, 2/40.0, 1/40.0 },
+	{ -1/40.0, -2/40.0, 0, 2/40.0, 1/40.0 },
+	{ -1/40.0, -2/40.0, 0, 2/40.0, 1/40.0 }
+};
+
+static const double filter_trivial[1][3] = {
+	{ -0.5, 0, 0.5 }
+};
+
 int main(int argc, char **argv)
 {
-	const char *infile, *outfile;
+	const char *infile, *outfile, *reffile;
 	double scale, offset;
-	int nmaplen;
-	unsigned char *nmapdata, *nmap;
+	int nmaplen, w, h;
+	unsigned char *nmapdata, *nmap, *refmap;
+	const char *filtertype;
+	const double *filter = NULL;
+	int filterw = 0, filterh = 0;
+#define USE_FILTER(f) \
+	do \
+	{ \
+		filterw = sizeof(*(f)) / sizeof(**(f)); \
+		filterh = sizeof((f)) / sizeof(*(f)); \
+		filter = &(f)[0][0]; \
+	} \
+	while(0)
 
 	if(argc > 1)
 		infile = argv[1];
@@ -685,14 +1042,24 @@ int main(int argc, char **argv)
 		return usage(*argv);
 	
 	if(argc > 3)
-		scale = atof(argv[3]);
+		filtertype = argv[3];
+	else
+		return usage(*argv);
+	
+	if(argc > 4)
+		scale = atof(argv[4]);
 	else
 		scale = 0;
 
-	if(argc > 4)
-		offset = atof(argv[4]);
+	if(argc > 5)
+		offset = atof(argv[5]);
+	else
+		offset = (scale<0) ? 1 : 0;
+
+	if(argc > 6)
+		reffile = argv[6];
 	else
-		offset = 0;
+		reffile = NULL;
 
 	nmapdata = FS_LoadFile(infile, &nmaplen);
 	if(!nmapdata)
@@ -707,7 +1074,56 @@ int main(int argc, char **argv)
 		printf("LoadTGA_BGRA failed\n");
 		return 2;
 	}
-	nmap_to_hmap(nmap, image_width, image_height, scale, offset);
+	w = image_width;
+	h = image_height;
+
+	if(reffile)
+	{
+		nmapdata = FS_LoadFile(reffile, &nmaplen);
+		if(!nmapdata)
+		{
+			printf("FS_LoadFile failed\n");
+			return 2;
+		}
+		refmap = LoadTGA_BGRA(nmapdata, nmaplen);
+		free(nmapdata);
+		if(!refmap)
+		{
+			printf("LoadTGA_BGRA failed\n");
+			return 2;
+		}
+		if(image_width != w || image_height != h)
+		{
+			printf("reference map must have same size as input normalmap\n");
+			return 2;
+		}
+	}
+	else
+		refmap = NULL;
+
+	if(!strcmp(filtertype, "trivial"))
+		USE_FILTER(filter_trivial);
+	if(!strcmp(filtertype, "prewitt3"))
+		USE_FILTER(filter_prewitt3);
+	if(!strcmp(filtertype, "scharr3"))
+		USE_FILTER(filter_scharr3);
+	if(!strcmp(filtertype, "sobel3"))
+		USE_FILTER(filter_sobel3);
+	if(!strcmp(filtertype, "prewitt5"))
+		USE_FILTER(filter_prewitt5);
+	if(!strcmp(filtertype, "sobel5"))
+		USE_FILTER(filter_sobel5);
+
+	if(scale < 0)
+	{
+		if(filter)
+			hmap_to_nmap_local(nmap, image_width, image_height, -scale-1, offset, filter, filterw, filterh);
+		else
+			hmap_to_nmap(nmap, image_width, image_height, -scale-1, offset);
+	}
+	else
+		nmap_to_hmap(nmap, refmap, image_width, image_height, scale, offset, filter, filterw, filterh);
+
 	if(!Image_WriteTGABGRA(outfile, image_width, image_height, nmap))
 	{
 		printf("Image_WriteTGABGRA failed\n");