diff options
Diffstat (limited to 'src/renderer.c')
-rw-r--r-- | src/renderer.c | 205 |
1 files changed, 205 insertions, 0 deletions
diff --git a/src/renderer.c b/src/renderer.c new file mode 100644 index 0000000..9384e30 --- /dev/null +++ b/src/renderer.c @@ -0,0 +1,205 @@ +/* + * renderer.c + * + * Render Fourier Transform output array into display array + * + * (c) 2006-2007 Thomas White <taw27@cam.ac.uk> + * + * synth2D - Two-Dimensional Crystallographic Fourier Synthesis + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdlib.h> +#include <string.h> +#include <stdio.h> +#include <math.h> +#include <fftw3.h> + +#include "renderer.h" + +double renderer_width(int width, int height, double gamma, int nx, int ny) { + return nx*width + ny*height*fabs(cos(gamma)); +} + +double renderer_height(int width, int height, double gamma, int nx, int ny) { + return height*sin(gamma)*ny; +} + +/* Provide the coordinate in the image (ie xn,yn) given x,y in the Fourier array */ +double renderer_map_x(double x, double y, int width, int height, double gamma, int nx, int ny) { + + double xn, offs; + + if ( gamma > M_PI_2 ) { + offs = (double)ny*height*cos(M_PI-gamma); + } else { + offs = 0.0; + } + xn = x + offs + y*cos(gamma); + + return xn; + +} + +/* Same, but return y */ +double renderer_map_y(double x, double y, int width, int height, double gamma, int nx, int ny) { + return y*sin(gamma); +} + +static double renderer_map_reverse_x(double xn, double yn, int width, int height, double gamma, int nx, int ny) { + + double offs, xd, yd; + + if ( gamma > M_PI_2 ) { + offs = (double)ny*height*cos(M_PI-gamma); + } else { + offs = 0.0; + } + + yd = (double)yn / sin(gamma); + xd = (double)xn - offs - yd*cos(gamma); + + while ( xd < 0 ) xd += width; + xd = fmod(xd, width); + + return xd; + +} + +static double renderer_map_reverse_y(double xn, double yn, int width, int height, double gamma, int nx, int ny) { + + double yd; + + yd = (double)yn / sin(gamma); + while ( yd < 0 ) yd += height; + yd = fmod(yd, height); + + return yd; + +} + +static double renderer_interpolate_linear_re(double xd, double yd, fftw_complex *out, int width, int height) { + + double frac, re, re1, re2; + int x, y; + + /* Get the left-hand point value */ + x = floor(xd); y = floor(yd); + if ( x >= width ) x -= width; + if ( y >= height ) y -= height; + re1 = out[(height-1-y)+height*(width-1-x)][0]; + + /* Get the right-hand point value */ + x++; + if ( x >= width ) x -= width; + re2 = out[(height-1-y)+height*(width-1-x)][0]; + + frac = fmod(xd, 1); + re = (1-frac)*re1 + frac*re2; + + return re; + +} + +static double renderer_interpolate_linear_im(double xd, double yd, fftw_complex *out, int width, int height) { + + double frac, im, im1, im2; + int x, y; + + /* Get the left-hand point value */ + x = floor(xd); y = floor(yd); + if ( x >= width ) x -= width; + if ( y >= height ) y -= height; + im1 = out[(height-1-y)+height*(width-1-x)][1]; + + /* Get the right-hand point value */ + x++; + if ( x >= width ) x -= width; + im2 = out[(height-1-y)+height*(width-1-x)][1]; + + frac = fmod(xd, 1); + im = (1-frac)*im1 + frac*im2; + + return im; + +} + +static double renderer_interpolate_bilinear_re(double xd, double yd, fftw_complex *out, int width, int height) { + + double frac, re, re1, re2; + + /* Get the lower interpolated value */ + re1 = renderer_interpolate_linear_re(xd, yd, out, width, height); + + /* Get the upper interpolated value */ + yd++; + if ( yd >= height ) yd -= height; + re2 = renderer_interpolate_linear_re(xd, yd, out, width, height); + + frac = fmod(yd, 1); + re = (1-frac)*re1 + frac*re2; + + return re; + +} + +static double renderer_interpolate_bilinear_im(double xd, double yd, fftw_complex *out, int width, int height) { + + double frac, im, im1, im2; + + /* Get the lower interpolated value */ + im1 = renderer_interpolate_linear_im(xd, yd, out, width, height); + + /* Get the upper interpolated value */ + yd++; + if ( yd >= height ) yd -= height; + im2 = renderer_interpolate_linear_im(xd, yd, out, width, height); + + frac = fmod(yd, 1); + im = (1-frac)*im1 + frac*im2; + + return im; + +} + +ComplexArray renderer_draw(fftw_complex *out, int width, int height, double gamma, int nx, int ny) { + + ComplexArray cxar; + int width_n, height_n; + int xn, yn; + + width_n = (int)renderer_width(width, height, gamma, nx, ny); + height_n = (int)renderer_height(width, height, gamma, nx, ny); + + cxar.re = malloc(height_n*width_n*sizeof(double)); + cxar.im = malloc(height_n*width_n*sizeof(double)); + + for ( yn=0; yn<height_n; yn++ ) { + for ( xn=0; xn<width_n; xn++ ) { + + double xd, yd; + double re, im; + + /* Map this pixel onto the FFTW output array */ + xd = renderer_map_reverse_x(xn, yn, width, height, gamma, nx, ny); + yd = renderer_map_reverse_y(xn, yn, width, height, gamma, nx, ny); + + /* Extract (interpolatively) the value */ + re = renderer_interpolate_bilinear_re(xd, yd, out, width, height); + im = renderer_interpolate_bilinear_im(xd, yd, out, width, height); + + /* Store the values */ + cxar.re[xn+width_n*yn] = re; + cxar.im[xn+width_n*yn] = im; + + } + } + + return cxar; + +} + |