#include #include #include #include #include #include /* Number of fractal rules. */ #define RULE_COUNT 3 /* Number of fractal iterations. */ #define ITERATIONS 16 /* Number of sub-iterations to apply in one step. Normal fractal works like this: 1. Select a rule at random, out of rule_count choices. 2. Apply rule. 3. Repeat steps 1-2 multiple times. We can speed this up like this: 0. Pre-compute explosion of rules_count**substep rules, i.e. all possible rules up to substep count. 1. Select one rule out of rule_count**substep choices. 2. Apply rule. 3. Repeat steps 1-2 multiple times. This is faster because the number of iterations at step 3 only need to be ITERATION_COUNT / SUBSTEP_COUNT, so we reduce the number of transformations per input pixel. */ #define SUBSTEP_COUNT 8 /* Maximum amount of overlap on the same pixel. */ #define MAX_OVERLAP 80 static png_image image; static png_bytep buffer[2]; static int s, i, j, k, w, h, out_of_bounds, *pixel, *source, /* Count of each pixel value. */ color_count[MAX_OVERLAP + 2], /* Output rule index. */ t; static double fx, fy, nx, ny, zx, zy, /* Base vertices for rules. The transforms will be centered around these vertices. */ vertex[RULE_COUNT][2], /* Transformation rules. Basically top 6 elements of a transformation matrix. [ x[0] x[1] x[2] x[3] x[4] x[5] 0.0 0.0 1.0 ] */ base_rule[RULE_COUNT][6], rule[6561][6], /* 6561 = 3 ** 8 */ identity[6] = {1, 0, 0, 0, 1, 0}; static double rnd(double range) { return (double)rand() / (RAND_MAX + 1u) * range; } /* Apply transformation a to b, storing result in b. */ static void ApplyTransformation(double *a, double *b) { double r[] = { a[0] * b[0] + a[1] * b[3], a[0] * b[1] + a[1] * b[4], a[0] * b[2] + a[1] * b[5] + a[2], a[3] * b[0] + a[4] * b[3], a[3] * b[1] + a[4] * b[4], a[3] * b[2] + a[4] * b[5] + a[5] }; memcpy(b, r, sizeof(r)); } /* Recursively compute all transformations. */ static void GenerateTransformations(int depth, double *parent) { double child[6]; int r; if( depth == SUBSTEP_COUNT ) memcpy(rule[t++], parent, 6 * sizeof(double)); else for(r = 0; r < RULE_COUNT; GenerateTransformations(depth + 1, child)) ApplyTransformation(base_rule[r++], (double*)memcpy(child, parent, 6 * sizeof(double))); } /* Compute area of parallelogram from 3 corner vertices. */ static double GetParallelogramArea(double ax, double ay, double bx, double by, double cx, double cy) { return fabs((bx - ax) * (cy - ay) - (by - ay) * (cx - ax)); } int main(int argc, char **argv) { if( argc != 4 ) return printf("%s \n", *argv); image.version = PNG_IMAGE_VERSION; if( png_image_begin_read_from_file(&image, argv[1]) ) { image.format = PNG_FORMAT_RGBA; for(i = 0; i < 2;) buffer[i++] = (png_bytep)calloc(s = (w = image.width) * (h = image.height), 4); k = *buffer && buffer[1] && png_image_finish_read(&image, NULL, *buffer, 0, NULL); } if( !k ) return printf("Error reading %s\n", argv[1]); srand(time(NULL)); /* Initialize vertices in [0.1,0.9] range with sufficient spread. */ for(zx = 0; zx < 0.4; zx = GetParallelogramArea(vertex[0][0], vertex[0][1], vertex[1][0], vertex[1][1], vertex[2][0], vertex[2][1])) for(i = 0; i < RULE_COUNT; i++) for(j = 0; j < 2;) vertex[i][j++] = rnd(0.8) + 0.1; /* Pick random points near the generated vertices. These will be used to initialize transformation rules. To ensure that each of these transformations produce interesting results, we also require that the parallelograms formed by the vertices cover sufficient area. */ for(i = 0; i < RULE_COUNT; i++) { if( sizeof('c') > 1 ) { /* C mode: generate constrained axis-aligned rules. */ base_rule[i][0] = fx = rnd(0.6) + 0.1 * (i + 2); base_rule[i][4] = fy = rnd(0.6) + 0.1 * (i + 2); base_rule[i][1] = base_rule[i][3] = 0.0; base_rule[i][2] = vertex[i][0] - fx * 0.5; base_rule[i][5] = vertex[i][1] - fy * 0.5; } else { /* C++ mode: generate free range random transforms that falls within expected area range. */ for(zx = 0; zx < 0.3 || zx > 0.7; zx = GetParallelogramArea(base_rule[i][0], base_rule[i][1], base_rule[i][2], base_rule[i][3], base_rule[i][4], base_rule[i][5])) for(j = 0; j < 6; j++) base_rule[i][j] = vertex[i][j % 2] + rnd(1) - 0.5; /* Convert vertices to transformation matrix. Assume that input specifies 3 vertices of a triangle, we want a transformation matrix such that F(0,0) -> A F(1,0) -> B F(0,1) -> C Transformations are applied as: X(x,y) = output->x[0] * x + output->x[1] * y + output->x[2] Y(x,y) = output->x[3] * x + output->x[4] * y + output->x[5] */ nx = base_rule[i][2]; ny = base_rule[i][3]; zx = base_rule[i][4]; zy = base_rule[i][5]; /* (0,0) -> A */ base_rule[i][2] = fx = base_rule[i][0]; base_rule[i][5] = fy = base_rule[i][1]; /* (1,0) -> B */ base_rule[i][0] = nx - fx; base_rule[i][3] = ny - fy; /* (0,1) -> C */ base_rule[i][1] = zx - fx; base_rule[i][4] = zy - fy; } } /* Pre-compute all transformations. */ GenerateTransformations(0, identity); /* Draw points. This is done by applying random rules iteratively to each starting point, and repeat the process until we appear to be getting diminishing returns, which is one of these conditions: - Many points are landing on pixels that were previously drawn. - Many points are drawn out of bounds. */ pixel = (int*)buffer[1]; for(out_of_bounds = 0; out_of_bounds < s * MAX_OVERLAP;) { fx = rnd(1); fy = rnd(1); for(i = 0; i < ITERATIONS / SUBSTEP_COUNT; i++) { k = (int)rnd(6561); nx = fx * rule[k][0] + fy * rule[k][1] + rule[k][2]; ny = fx * rule[k][3] + fy * rule[k][4] + rule[k][5]; fx = nx; fy = ny; } /* Convert fractal coordinates from [0,1] range back to screen range. */ i = (int)round(fx * w); j = (int)round(fy * h); if( i < 0 || j < 0 || i >= w || j >= h ) out_of_bounds++; else if( pixel[j * w + i]++ >= MAX_OVERLAP ) break; /* Too many pixels landed on the same point. */ } /* Count number of colors, and pick a threshold such that number of pixels above and below the threshold are roughly equal. */ for(i = j = 0; i < s;) color_count[pixel[i++]]++; for(k = MAX_OVERLAP; k > 1 && (j += color_count[k]) < s / 2; k--); /* Move pixels according to threshold. */ source = (int*)*buffer; for(i = 0; i++ < s; pixel++, *source++ *= j) *pixel = (j = *pixel < k) ? 0 : *source; for(i = 0; i < 2; i++) if( !png_image_write_to_file(&image, argv[i + 2], 0, buffer[i], 0, NULL) ) return printf("Error writing %s\n", argv[i + 2]); return 0; }