void nmap_to_hmap(unsigned char *map, const unsigned char *refmap, int w, int h, double scale, double offset)
{
int x, y;
+ int fx, fy;
+ double ffx, ffy;
double nx, ny, nz;
double v, vmin, vmax;
#ifndef C99
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)
for(y = 0; y < h; ++y)
for(x = 0; x < w; ++x)
{
- int fx = x;
- int fy = y;
+ fx = x;
+ fy = y;
if(fx > w/2)
fx -= w;
if(fy > h/2)
fy -= h;
+ /* 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 * (fx * freqspace1[(w*y+x)] + fy * freqspace2[(w*y+x)]) / (fx*fx + fy*fy) / TWO_PI;
+ 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] = -(fx * freqspace1[(w*y+x)][1] + fy * freqspace2[(w*y+x)][1]) / (fx*fx + fy*fy) / TWO_PI;
- freqspace1[(w*y+x)][1] = (fx * save + fy * freqspace2[(w*y+x)][0]) / (fx*fx + fy*fy) / TWO_PI;
+ 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
{
fftw_execute(f12i1);
- /* renormalize */
+ /* renormalize, find min/max */
+ vmin = vmax = 0;
for(y = 0; y < h; ++y)
for(x = 0; x < w; ++x)
{
#ifdef C99
- imgspace1[(w*y+x)] /= (w*h);
+ v = creal(imgspace1[(w*y+x)] /= (w*h));
#else
- imgspace1[(w*y+x)][0] /= (w*h);
+ 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)
s = 1;
}
- printf("Ref-computed scale: %f\nRef-computed offset: %f\n", s, o);
-
/*
* now apply user-given offset and scale to these values
* (x * s + o) * scale + offset
}
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;
- }
-
/*
* 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);
}
+ 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)