fern.txt

64
1024
4
0.20 -0.26 0.23 0.22 0 1.60 0.07 
-0.15 0.28 0.26 0.24 0 0.44 0.07 0.85 
0.04 -0.04 0.85 0 1.60 0.85 
0 0 0 0.16 0 0 0.01 

ifs.cc

#include <fstream>
#include <iostream>
#include <cstdlib>
#include <complex>

int main(int argc, char* argv[]) {
	std::ifstream infile;
	infile.open(argv[1], std::ifstream::in);
	int N;
	const int N1 = 1 << 18;
	int m;
	int n;

	infile >> m;
	infile >> n;
	infile >> N;

	std::cerr << "m : " << m << std::endl;
	std::cerr << "n : " << n << std::endl;
	std::cerr << "N : " << N << std::endl;

	double** a;
	double* p;
	a =  new double*[N];
	p = new double[N];

	for(int i = 0; i < N; i++) {
		a[i] = new double[6];

		for(int j = 0; j < 6; j++) {
			double t;
			infile >> t;
			a[i][j] = t;
		}

		infile >> p[i];

	}
	infile.close();

	for(int i = 1; i < N; i++) {
		p[i] += p[i - 1];
	}

	for(int i = 0; i < N; i++) {
		for(int j = 0; j < 2; j++) {
			std::cerr << "[ ";

			for(int k = 0; k < 2; k++) {
				std::cerr << a[i][j * 2 + k];
				std::cerr << " ";
			}

			std::cerr << "] [ ";
			std::cerr << a[i][j + 4];
			std::cerr << " ]" << std::endl;
		}

		std::cerr << std::endl;
	}

	double x, y;
	x = 0;
	y = 0;
	srand48(time(NULL));
	double left = 0, right = 0, top = 0, bottom = 0;

	for(long long i = 0; i < N1; i++) {
		double x2, y2;

		double prob = drand48();
		int idx;

		for(idx = 0; idx < N - 1; idx++) {
			if(prob < p[idx]) {
				break;
			}
		}

		x2 = a[idx][0] * x + a[idx][1] * y + a[idx][4];
		y2 = a[idx][2] * x + a[idx][3] * y + a[idx][5];
		x = x2;
		y = y2;

		left = std::min(x, left);
		bottom = std::min(y, bottom);
		right = std::max(x, right);
		top = std::max(y, top);
	}

	const double w = 0.02 * (((right - left) + (top - bottom)) / 2);
	left -= w;
	bottom -= w;
	top += w;
	right += w;


	double width = right - left;
	double height = top - bottom;
	long long n0, n1;

	if(width > height) {
		n0 = width * n / height;
		n1 = n;
	} else {
		n0 = n;
		n1 = height * n / width;
	}

	std::cerr << n0 << "," << n1 << ":" << m << std::endl;
	std::cerr << left << "," << bottom << std::endl;
	std::cerr << right << "," << top << std::endl;

	int* counts = new int[n0 * n1];
	long long iterations = n0 * n1 * m;

	for(long long i = 0; i < iterations; i++) {
		double x2, y2;
		double prob = drand48();
		int idx;

		for(idx = 0; idx < N - 1; idx++) {
			if(prob < p[idx]) {
				break;
			}
		}

		x2 = a[idx][0] * x + a[idx][1] * y + a[idx][4];
		y2 = a[idx][2] * x + a[idx][3] * y + a[idx][5];
		x = x2;
		y = y2;
		int x_s = n0 * (x - left) / width;
		int y_s = n1 * (top - y) / height;

		if(0 <= x_s && x_s < n0 && 0 <= y_s && y_s < n1) {
			counts[x_s + y_s * n0] ++;
		}
	}

	std::ofstream outfile;
	outfile.open(argv[2]);
	{
		int t0 = sizeof(counts[0]);
		outfile.write(reinterpret_cast<const char*>(&t0), sizeof(int));
	}
	outfile.write(reinterpret_cast<const char*>(&n0), sizeof(int));
	outfile.write(reinterpret_cast<const char*>(&n1), sizeof(int));
	outfile.write(reinterpret_cast<const char*>(counts), sizeof(int)*n0 * n1);
	outfile.close();

	delete [] counts;

	for(int i = 0; i < N; i++) {
		delete [] a[i];
	}

	delete [] a;
	delete [] p;
	return 0;
}

maple.txt

64
2048
4
0.14	0.01	0.00	0.51	-0.08	-1.31	0.10
0.43	0.52	-0.45	0.50	1.49	-0.75	0.35
0.45	-0.49	0.47	0.47	-1.62	-0.74	0.35
0.49	0.00	0.00	0.51	0.02	1.62	0.20

process_ifs.py

def main():
    import numpy as np
    import sys
    in_fname = sys.argv[1]
    out_fname = sys.argv[2]
    with open(in_fname, 'rb') as infile:
        t = int.from_bytes(infile.read(4), byteorder='little')
        n0 = int.from_bytes(infile.read(4), byteorder='little')
        n1 = int.from_bytes(infile.read(4), byteorder='little')
        if t == 4:
            arr = np.fromfile(infile, dtype=np.int32)
        elif t == 8:
            arr = np.fromfile(infile, dtype=np.int64)
        else:
            raise
        print(n0, n1)
    img = arr.reshape((n1, n0)).astype(np.float64)
    img -= np.min(img)
    img /= np.max(img)

    img = pow(img, 1 / 7)

    img -= np.min(img)
    img /= np.max(img)
    img = 1 - img
    img = (img * 255).astype(np.uint8)
    import imageio
    imageio.imsave(out_fname, img)


if __name__ == '__main__':
    main()