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)
fy -= h;
#ifdef C99
if(fx||fy)
- freqspace1[(w*y+x)] = I * (fx * freqspace1[(w*y+x)] + fy * freqspace2[(w*y+x)]) / (fx*fx + fy*fy) / TWO_PI;
+ freqspace1[(w*y+x)] = _Complex_I * (fx * freqspace1[(w*y+x)] + fy * freqspace2[(w*y+x)]) / (fx*fx + fy*fy) / TWO_PI;
else
freqspace1[(w*y+x)] = 0;
#else
fftw_execute(f12i1);
- if(refmap)
+ /* renormalize, find min/max */
+ vmin = vmax = 0;
+ for(y = 0; y < h; ++y)
+ for(x = 0; x < w; ++x)
{
- // refmap: a reference map to define the heights
- // alpha = weight, color = value
- // if more than one color value is used, colors are also matched
-
- // we do linear regression, basically
- // f'(x, y) = f(x, y) * scale + offset
- // sum((f(x, y) * scale + offset - ref_y(x, y))^2 * ref_a(x, y)) minimize
-
- // diff by offset:
- // sum(-2*ref_y(x,y)*ref_a(x,y) + 2*scale*f(x,y)*ref_a(x,y) + 2*offset*ref_a(x,y)) = 0
- // diff by scale:
- // sum(-2*f(x,y)*ref_a(x,y) + 2*scale*f(x,y)^2*ref_a(x,y) + 2*offset*f(x,y)*ref_a(x,y)) = 0
- // ->
- // offset = (sfa*sfya - sffa*sya) / (sfa*sfa-sa*sffa)
- // scale = (sfa*sya - sa*sfya) / (sfa*sfa-sa*sffa)
+#ifdef C99
+ v = creal(imgspace1[(w*y+x)] /= (w*h));
+#else
+ v = (imgspace1[(w*y+x)][0] /= (w*h));
+ imgspace1[(w*y+x)][1] /= (w*h);
+#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;
for(y = 0; y < h; ++y)
for(x = 0; x < w; ++x)
{
- a = (int)refmap[(w*y+x)*4+0];
- 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);
- v = (v - 128.0) / 127.0; // value 0 is forbidden, 1 -> -1, 255 -> 1
+ 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
f = creal(imgspace1[(w*y+x)]);
#else
#endif
if(a <= 0)
continue;
- if(y < mi)
- mi = y;
- if(y > ma)
- ma = y;
+ 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;
}
- sfa /= (w*h);
- sffa /= (w*h);
- sffa /= (w*h);
- sfva /= (w*h);
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
+ else /* all values of v are equal, so we cannot get scale; we can still get offset */
{
o = ((sva - sfa) / sa);
s = 1;
}
- // now apply user-given offset and scale to these values
- // (x * s + o) * scale + offset
- // x * s * scale + o * scale + offset
+
+ /*
+ * 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)
{
-#ifdef C99
- vmin = vmax = creal(imgspace1[0]);
-#else
- vmin = vmax = imgspace1[0][0];
-#endif
- for(y = 0; y < h; ++y)
- for(x = 0; x < w; ++x)
- {
-#ifdef C99
- v = creal(imgspace1[(w*y+x)]);
-#else
- v = imgspace1[(w*y+x)][0];
-#endif
- if(v < vmin)
- vmin = v;
- if(v > vmax)
- vmax = v;
- }
-
- vmin /= (w*h);
- vmax /= (w*h);
-
/*
* 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)
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 f12i1 = fftw_plan_dft_2d(w, h, freqspace1, imgspace1, FFTW_BACKWARD, FFTW_ESTIMATE);
- fftw_plan f22i2 = fftw_plan_dft_2d(w, h, freqspace2, imgspace2, FFTW_BACKWARD, FFTW_ESTIMATE);
+ 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)
freqspace1[(w*y+x)] *= 1 - pow(abs(fx) / (double)(w/2), 1);
freqspace1[(w*y+x)] *= 1 - pow(abs(fy) / (double)(h/2), 1);
- freqspace2[(w*y+x)] = TWO_PI*I * fy * freqspace1[(w*y+x)]; /* y derivative */
- freqspace1[(w*y+x)] = TWO_PI*I * fx * freqspace1[(w*y+x)]; /* x derivative */
+ 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
/* a lowpass to prevent the worst */
freqspace1[(w*y+x)][0] *= 1 - pow(abs(fx) / (double)(w/2), 1);
if(reffile)
{
- nmapdata = FS_LoadFile(infile, &nmaplen);
+ nmapdata = FS_LoadFile(reffile, &nmaplen);
if(!nmapdata)
{
printf("FS_LoadFile failed\n");
return 2;
}
}
+ else
+ refmap = NULL;
if(scale < -6)
hmap_to_nmap(nmap, image_width, image_height, -scale-7, offset);