More filters

1 files changed, 133 insertions, 7 deletions

diff --git a/jinc.c b/jinc.c
index 3eff78c..a0f53d3 100644
--- a/jinc.c
+++ b/jinc.c
@@ -23,7 +23,10 @@ static const double J1_roots[] =
 static double
 jinc (double d)
 {
-    return j1 (M_PI * d) / (2 * d);
+    if (d == 0.0)
+        return AT_ZERO;
+    else
+        return j1 (M_PI * d) / (2 * d);
 }
 
 /* Jinc windowed jinc function */
@@ -33,9 +36,7 @@ jinc_jinc (double x, double y)
     double d = sqrt (x * x + y * y);
     double rf;
 
-    if (d == 0.0)
-        rf = AT_ZERO;
-    else if (fabs (d) <= LOBES)
+    if (fabs (d) <= LOBES)
         rf = jinc (d) * jinc (d * (MAIN_LOBE / LOBES));
     else
         rf = 0;
@@ -43,7 +44,7 @@ jinc_jinc (double x, double y)
     return rf;
 }
 
-#define SIZE 256
+#define SIZE 512
 
 static complex_image_t *
 make_jinc_image (void)
@@ -63,7 +64,7 @@ make_jinc_image (void)
 	    double y = - LOBES + step / 2 + j * step;
 	    double v;
 
-	    v = jinc_jinc (x, y);
+	    v = jinc_jinc (8 * x, 8 * y);
 
 	    v = 0.5 * (v + 1);
 	    
@@ -78,6 +79,131 @@ make_jinc_image (void)
     return result;
 }
 
+static double
+sinc (double x)
+{
+    return sin (M_PI * x) / (M_PI * x);
+}
+
+#define LLOBES (3)
+
+static double
+lanczos (double x)
+{
+    if (fabs (x) < LLOBES)
+        return sinc (x) * sinc (x / LLOBES);
+    else
+        return 0.0;
+}
+
+static double
+lanczos2d (double x, double y)
+{
+    return lanczos (x) * lanczos (y);
+}
+
+static complex_image_t *
+make_lanczos_image (void)
+{
+    complex_image_t *result = complex_image_new (SIZE, SIZE);
+    int i, j;
+
+    for (j = 0; j < SIZE; ++j)
+    {
+	for (i = 0; i < SIZE; ++i)
+	{
+	    complex_t *r = &(result->red[j * SIZE + i]);
+	    complex_t *g = &(result->green[j * SIZE + i]);
+	    complex_t *b = &(result->blue[j * SIZE + i]);
+	    double step = (2 * LLOBES) / (double)SIZE;
+	    double x = - LLOBES + step / 2 + i * step;
+	    double y = - LLOBES + step / 2 + j * step;
+	    double v;
+
+	    v = lanczos2d (8 * x, 8 * y);
+
+	    v = 0.5 * (v + 1);
+	    
+	    r->re = g->re = b->re = 0.0;
+	    r->re = 1 - v;
+            g->re = v;
+            b->re = 0;
+
+	    r->im = g->im = b->im = 0.0;
+	}
+    }
+
+    return result;
+}
+
+static complex_image_t *
+make_sinc_image (void)
+{
+    complex_image_t *result = complex_image_new (SIZE, SIZE);
+    int i, j;
+
+    for (j = 0; j < SIZE; ++j)
+    {
+	for (i = 0; i < SIZE; ++i)
+	{
+	    complex_t *r = &(result->red[j * SIZE + i]);
+	    complex_t *g = &(result->green[j * SIZE + i]);
+	    complex_t *b = &(result->blue[j * SIZE + i]);
+	    double step = (2 * LLOBES) / (double)SIZE;
+	    double x = - LLOBES + step / 2 + i * step;
+	    double y = - LLOBES + step / 2 + j * step;
+	    double v;
+
+	    v = sinc (8 * x) * sinc (8 * y);
+
+	    v = 0.5 * (v + 1);
+	    
+	    r->re = g->re = b->re = 0.0;
+	    r->re = 1 - v;
+            g->re = v;
+            b->re = 0;
+
+	    r->im = g->im = b->im = 0.0;
+	}
+    }
+
+    return result;
+}
+
+static complex_image_t *
+make_jnw_image (void)
+{
+    complex_image_t *result = complex_image_new (SIZE, SIZE);
+    int i, j;
+
+    for (j = 0; j < SIZE; ++j)
+    {
+	for (i = 0; i < SIZE; ++i)
+	{
+	    complex_t *r = &(result->red[j * SIZE + i]);
+	    complex_t *g = &(result->green[j * SIZE + i]);
+	    complex_t *b = &(result->blue[j * SIZE + i]);
+	    double step = (2 * LOBES) / (double)SIZE;
+	    double x = - LOBES + step / 2 + i * step;
+	    double y = - LOBES + step / 2 + j * step;
+	    double v;
+
+	    v = jinc (sqrt (8 * x * 8 * x + 8 * y * 8 * y));
+
+	    v = 0.5 * (v + 1);
+	    
+	    r->re = g->re = b->re = 0.0;
+	    b->re = v;
+	    g->re = 1 - v;
+            r->re = v;
+
+	    r->im = g->im = b->im = 0.0;
+	}
+    }
+
+    return result;
+}
+
 static complex_t
 filter (complex_t c, double dist)
 {
@@ -120,7 +246,7 @@ main (int argc, char **argv)
     
     g_type_init ();
     
-    image = make_jinc_image ();
+    image = make_jnw_image ();
     
     g_print ("width, height %d %d\n", image->width, image->height);