aboutsummaryrefslogtreecommitdiff
path: root/src/renderer.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/renderer.c')
-rw-r--r--src/renderer.c205
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;
+
+}
+