use pgm

Author: mattip

create_fractal.c

@@ -1,6 +1,6 @@
-int mandel(int x, int y, int max_iters, unsigned char * val);
+int mandel(int x, int y, int max_iters, unsigned int * val);
 
-void create_fractal(unsigned char * image, int width, int height, 
+void create_fractal(unsigned int * image, int width, int height, 
                     int iters) {
     float pixel_size_x = 3.0 / width;
     float pixel_size_y = 2.0 / height;
@@ -9,7 +9,7 @@ void create_fractal(unsigned char * image, int width, int height,
         int yy = y * width;
         for (int x=0; x<width; x++) {
             float real = x * pixel_size_x - 2;
-            unsigned char color;
+            unsigned int color;
             int ret = mandel(real, imag, iters, &color);
             image[yy + x] = color;
         }

create_fractal.py

@@ -0,0 +1,28 @@
+from __future__ import print_function, division
+
+from mandel import mandel
+
+def create_fractal(image, width, height, iters):
+    pixel_size_x = 3.0 / width
+    pixel_size_y = 2.0 / height
+    for y in range(height):
+        imag = y * pixel_size_y - 1
+        yy = y*width
+        for x in range(width):
+            real = x * pixel_size_x - 2
+            color = mandel(real, imag, iters)
+            image[yy+x] = color
+
+if __name__ == '__main__':
+    from timeit import default_timer as timer
+    import numpy as np
+    from PIL import Image
+    width = 1500
+    height = 1000
+    image = np.empty(width*height, dtype=np.uint8)
+    s = timer()
+    create_fractal(image, width, height, 20)
+    e = timer()
+    print('pure python required {}f.2 secs'.format(e - s))
+    im = Image.frombuffer('L', (width, height), image, 'raw', 'L', 0, 1)
+    im.save('python_numpy.png') 

main.c

@@ -1,9 +1,25 @@
 #include <stdio.h>
-#include <pbm.h>
+#include <pgm.h>
+#include <time.h>
 
-int mandel(int x, int y, int max_iters, unsigned char * val);
+// call this function to start a nanosecond-resolution timer
+struct timespec timer_start(){
+    struct timespec start_time;
+    clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &start_time);
+    return start_time;
+}
+
+// call this function to end a timer, returning nanoseconds elapsed as a long
+long timer_end(struct timespec start_time){
+    struct timespec end_time;
+    clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &end_time);
+    long diffInNanos = end_time.tv_nsec - start_time.tv_nsec;
+    return diffInNanos;
+}
+
+int mandel(int x, int y, int max_iters, unsigned int * val);
 
-void create_fractal(unsigned char * image, int width, int height, 
+void create_fractal(unsigned int * image, int width, int height, 
                     int iters);
 
 int main(int argc, const char *argv[], const char * env[])
@@ -12,16 +28,21 @@ int main(int argc, const char *argv[], const char * env[])
     int height = 1000;
     int iters = 20;
     FILE * fid = NULL;
-    unsigned char * image = (unsigned char*)malloc(width*height);
+    struct timespec vartime;
+    long time_elapsed_nanos;
+    unsigned int * image = (unsigned int*)malloc(width * height * sizeof(unsigned int));
     if (NULL == image)
         return -1;
-    fid = fopen("c.pbm", "wb");
+    fid = fopen("c.pgm", "wb");
     if (NULL == fid)
         return -2;
+    vartime = timer_start();
     create_fractal(image, width, height, iters);
-    pbm_writepbminit(fid, width, height, 0);
+    time_elapsed_nanos = timer_end(vartime);
+    fprintf(stdout, "create_fractal required %ld millisecs\n", time_elapsed_nanos / 1000000);
+    pgm_writepgminit(fid, width, height, 255, 0);
     for (int i=0; i<height; i++)
-        pbm_writepbmrow(fid, image + width*i, width, 0);
+        pgm_writepgmrow(fid, image + width*i, width, 255, 0);
     fclose(fid);
     return 0;
 }

mandel.c

@@ -18,7 +18,7 @@ def mandel(x, y, max_iters):
 
 #include <complex.h>
 
-int mandel(int x, int y, int max_iters, unsigned char * val)
+int mandel(int x, int y, int max_iters, unsigned int * val)
 {
     int i = 0;
     _Complex float c = CMPLX(x, y);

mandel.py

@@ -1,57 +1,5 @@
-#! /usr/bin/env python
-# -*- coding: utf-8 -*-
-from __future__ import print_function, division, absolute_import
 
-'''
-Taken from https://github.com/numba/numba/blob/master/examples/mandel/mandel_jit.py
-On June 3, 2017, then modified by mattip
-Original copyright:
-
-Copyright (c) 2012, Continuum Analytics, Inc.
-All rights reserved.
-
-Redistribution and use in source and binary forms, with or without
-modification, are permitted provided that the following conditions are
-met:
-
-Redistributions of source code must retain the above copyright notice,
-this list of conditions and the following disclaimer.
-
-Redistributions in binary form must reproduce the above copyright
-notice, this list of conditions and the following disclaimer in the
-documentation and/or other materials provided with the distribution.
-THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
-LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
-A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
-HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
-SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
-LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
-DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
-THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
-(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
-OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-
-    Contact GitHub API Training Shop Blog About 
-
-'''
-
-from timeit import default_timer as timer
-
-import sys
-fname = 'test'
-if len(sys.argv) > 1 and 'nojit' in sys.argv[1]:
-   def jit(f):
-        print('pure python')
-        return f
-   fname = 'test_nojit'
-else:
-   from numba import jit
-   fname = 'test_numba'
-
-
-@jit
-def mandel(x, y, max_iters):
+def mandel(x, y, max_iters, value):
     """
     Given the real and imaginary parts of a complex number,
     determine if it is a candidate for membership in the Mandelbrot
@@ -63,51 +11,8 @@ def mandel(x, y, max_iters):
     for i in range(max_iters):
         z = z*z + c
         if (z.real*z.real + z.imag*z.imag) >= 4:
-            return i
-
-    return 255
-
-@jit
-def create_fractal(image, width, height, iters,
-                   min_x=-2.0, max_x=1.0, min_y=-1.0, max_y=1.0):
-    pixel_size_x = (max_x - min_x) / width
-    pixel_size_y = (max_y - min_y) / height
-    for y in range(height):
-        imag = min_y + y * pixel_size_y
-        yy = y*width
-        for x in range(width):
-            real = min_x + x * pixel_size_x
-            color = mandel(real, imag, iters)
-            image[yy+x] = color
-
-    return image
+            value[0] = i
+            return 0
+    value[i] = 255
+    return 1
 
-width = 1500
-height = 1000
-if 'cffi' in sys.argv:
-    import cffi
-    ffi = cffi.FFI()
-    image = ffi.new('unsigned char[{}]'.format(width*height))
-    img_as_bytes = bytes(ffi.buffer(image, width*height))
-    fname += '_cffi'
-else:
-    import numpy as np
-    image = np.empty(width*height, dtype=np.uint8)
-    fname += '_numpy'
-    img_as_bytes = image
-s = timer()
-create_fractal(image, width, height, 20)
-e = timer()
-print(e - s)
-if 0:
-    from matplotlib.pylab import imshow, jet, show, ion
-    imshow(image)
-    show()
-else:
-    from PIL import Image
-    try:
-        im = Image.frombuffer('L', (width, height), img_as_bytes, 'raw', 'L', 0, 1)
-        im.save(fname + '.png')
-    except:
-        import traceback;traceback.print_exc()
-        import pdb;pdb.set_trace()