|
| 1 | +# Copyright (c) 2020, Oracle and/or its affiliates. |
| 2 | +# Copyright (c) 2017, The PyPy Project |
| 3 | +# |
| 4 | +# The MIT License |
| 5 | +# Permission is hereby granted, free of charge, to any person |
| 6 | +# obtaining a copy of this software and associated documentation |
| 7 | +# files (the "Software"), to deal in the Software without |
| 8 | +# restriction, including without limitation the rights to use, |
| 9 | +# copy, modify, merge, publish, distribute, sublicense, and/or |
| 10 | +# sell copies of the Software, and to permit persons to whom the |
| 11 | +# Software is furnished to do so, subject to the following conditions: |
| 12 | +# |
| 13 | +# The above copyright notice and this permission notice shall be included |
| 14 | +# in all copies or substantial portions of the Software. |
| 15 | +# |
| 16 | +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS |
| 17 | +# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
| 18 | +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL |
| 19 | +# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
| 20 | +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING |
| 21 | +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER |
| 22 | +# DEALINGS IN THE SOFTWARE. |
| 23 | + |
| 24 | +# Copyright (C) 2005, Carl Friedrich Bolz |
| 25 | + |
| 26 | +"""create chaosgame-like fractals |
| 27 | +""" |
| 28 | + |
| 29 | +from __future__ import division, print_function |
| 30 | + |
| 31 | +import operator |
| 32 | +import random |
| 33 | +import math |
| 34 | +from functools import reduce |
| 35 | + |
| 36 | +random.seed(1234) |
| 37 | + |
| 38 | +class GVector(object): |
| 39 | + def __init__(self, x = 0, y = 0, z = 0): |
| 40 | + self.x = x |
| 41 | + self.y = y |
| 42 | + self.z = z |
| 43 | + |
| 44 | + def Mag(self): |
| 45 | + return math.sqrt(self.x ** 2 + self.y ** 2 + self.z ** 2) |
| 46 | + |
| 47 | + def dist(self, other): |
| 48 | + return math.sqrt((self.x - other.x) ** 2 + |
| 49 | + (self.y - other.y) ** 2 + |
| 50 | + (self.z - other.z) ** 2) |
| 51 | + |
| 52 | + def __add__(self, other): |
| 53 | + if not isinstance(other, GVector): |
| 54 | + raise ValueError( "Can't add GVector to " + str(type(other))) |
| 55 | + v = GVector(self.x + other.x, self.y + other.y, self.z + other.z) |
| 56 | + return v |
| 57 | + |
| 58 | + def __sub__(self, other): |
| 59 | + return self + other * -1 |
| 60 | + |
| 61 | + def __mul__(self, other): |
| 62 | + v = GVector(self.x * other, self.y * other, self.z * other) |
| 63 | + return v |
| 64 | + __rmul__ = __mul__ |
| 65 | + |
| 66 | + def linear_combination(self, other, l1, l2=None): |
| 67 | + if l2 is None: |
| 68 | + l2 = 1 - l1 |
| 69 | + v = GVector(self.x * l1 + other.x * l2, |
| 70 | + self.y * l1 + other.y * l2, |
| 71 | + self.z * l1 + other.z * l2) |
| 72 | + return v |
| 73 | + |
| 74 | + |
| 75 | + def __str__(self): |
| 76 | + return "<%f, %f, %f>" % (self.x, self.y, self.z) |
| 77 | + |
| 78 | + def __repr__(self): |
| 79 | + return "GVector(%f, %f, %f)" % (self.x, self.y, self.z) |
| 80 | + |
| 81 | +def GetKnots(points, degree): |
| 82 | + knots = [0] * degree + range(1, len(points) - degree) |
| 83 | + knots += [len(points) - degree] * degree |
| 84 | + return knots |
| 85 | + |
| 86 | +class Spline(object): |
| 87 | + """Class for representing B-Splines and NURBS of arbitrary degree""" |
| 88 | + def __init__(self, points, degree = 3, knots = None): |
| 89 | + """Creates a Spline. points is a list of GVector, degree is the |
| 90 | +degree of the Spline.""" |
| 91 | + if knots == None: |
| 92 | + self.knots = GetKnots(points, degree) |
| 93 | + else: |
| 94 | + if len(points) > len(knots) - degree + 1: |
| 95 | + raise ValueError("too many control points") |
| 96 | + elif len(points) < len(knots) - degree + 1: |
| 97 | + raise ValueError("not enough control points") |
| 98 | + last = knots[0] |
| 99 | + for cur in knots[1:]: |
| 100 | + if cur < last: |
| 101 | + raise ValueError("knots not strictly increasing") |
| 102 | + last = cur |
| 103 | + self.knots = knots |
| 104 | + self.points = points |
| 105 | + self.degree = degree |
| 106 | + |
| 107 | + def GetDomain(self): |
| 108 | + """Returns the domain of the B-Spline""" |
| 109 | + return (self.knots[self.degree - 1], |
| 110 | + self.knots[len(self.knots) - self.degree]) |
| 111 | + |
| 112 | + def __call__(self, u): |
| 113 | + """Calculates a point of the B-Spline using de Boors Algorithm""" |
| 114 | + dom = self.GetDomain() |
| 115 | + if u < dom[0] or u > dom[1]: |
| 116 | + raise ValueError("Function value not in domain") |
| 117 | + if u == dom[0]: |
| 118 | + return self.points[0] |
| 119 | + if u == dom[1]: |
| 120 | + return self.points[-1] |
| 121 | + I = self.GetIndex(u) |
| 122 | + d = [self.points[I - self.degree + 1 + ii] |
| 123 | + for ii in range(self.degree + 1)] |
| 124 | + U = self.knots |
| 125 | + for ik in range(1, self.degree + 1): |
| 126 | + for ii in range(I - self.degree + ik + 1, I + 2): |
| 127 | + ua = U[ii + self.degree - ik] |
| 128 | + ub = U[ii - 1] |
| 129 | + co1 = (ua - u) / (ua - ub) |
| 130 | + co2 = (u - ub) / (ua - ub) |
| 131 | + index = ii - I + self.degree - ik - 1 |
| 132 | + d[index] = d[index].linear_combination(d[index + 1], co1, co2) |
| 133 | + return d[0] |
| 134 | + |
| 135 | + def GetIndex(self, u): |
| 136 | + dom = self.GetDomain() |
| 137 | + for ii in range(self.degree - 1, len(self.knots) - self.degree): |
| 138 | + if u >= self.knots[ii] and u < self.knots[ii + 1]: |
| 139 | + I = ii |
| 140 | + break |
| 141 | + else: |
| 142 | + I = dom[1] - 1 |
| 143 | + return I |
| 144 | + |
| 145 | + def __len__(self): |
| 146 | + return len(self.points) |
| 147 | + |
| 148 | + def __repr__(self): |
| 149 | + return "Spline(%r, %r, %r)" % (self.points, self.degree, self.knots) |
| 150 | + |
| 151 | + |
| 152 | +class Chaosgame(object): |
| 153 | + def __init__(self, splines, thickness=0.1): |
| 154 | + self.splines = splines |
| 155 | + self.thickness = thickness |
| 156 | + self.minx = min([p.x for spl in splines for p in spl.points]) |
| 157 | + self.miny = min([p.y for spl in splines for p in spl.points]) |
| 158 | + self.maxx = max([p.x for spl in splines for p in spl.points]) |
| 159 | + self.maxy = max([p.y for spl in splines for p in spl.points]) |
| 160 | + self.height = self.maxy - self.miny |
| 161 | + self.width = self.maxx - self.minx |
| 162 | + self.num_trafos = [] |
| 163 | + maxlength = thickness * self.width / self.height |
| 164 | + for spl in splines: |
| 165 | + length = 0 |
| 166 | + curr = spl(0) |
| 167 | + for i in range(1, 1000): |
| 168 | + last = curr |
| 169 | + t = 1 / 999 * i |
| 170 | + curr = spl(t) |
| 171 | + length += curr.dist(last) |
| 172 | + self.num_trafos.append(max(1, int(length / maxlength * 1.5))) |
| 173 | + self.num_total = reduce(operator.add, self.num_trafos, 0) |
| 174 | + |
| 175 | + |
| 176 | + def get_random_trafo(self): |
| 177 | + r = random.randrange(int(self.num_total) + 1) |
| 178 | + l = 0 |
| 179 | + for i in range(len(self.num_trafos)): |
| 180 | + if r >= l and r < l + self.num_trafos[i]: |
| 181 | + return i, random.randrange(self.num_trafos[i]) |
| 182 | + l += self.num_trafos[i] |
| 183 | + return len(self.num_trafos) - 1, random.randrange(self.num_trafos[-1]) |
| 184 | + |
| 185 | + def transform_point(self, point, trafo=None): |
| 186 | + x = (point.x - self.minx) / self.width |
| 187 | + y = (point.y - self.miny) / self.height |
| 188 | + if trafo is None: |
| 189 | + trafo = self.get_random_trafo() |
| 190 | + start, end = self.splines[trafo[0]].GetDomain() |
| 191 | + length = end - start |
| 192 | + seg_length = length / self.num_trafos[trafo[0]] |
| 193 | + t = start + seg_length * trafo[1] + seg_length * x |
| 194 | + basepoint = self.splines[trafo[0]](t) |
| 195 | + if t + 1/50000 > end: |
| 196 | + neighbour = self.splines[trafo[0]](t - 1/50000) |
| 197 | + derivative = neighbour - basepoint |
| 198 | + else: |
| 199 | + neighbour = self.splines[trafo[0]](t + 1/50000) |
| 200 | + derivative = basepoint - neighbour |
| 201 | + if derivative.Mag() != 0: |
| 202 | + basepoint.x += derivative.y / derivative.Mag() * (y - 0.5) * \ |
| 203 | + self.thickness |
| 204 | + basepoint.y += -derivative.x / derivative.Mag() * (y - 0.5) * \ |
| 205 | + self.thickness |
| 206 | + else: |
| 207 | + print("r", end='') |
| 208 | + self.truncate(basepoint) |
| 209 | + return basepoint |
| 210 | + |
| 211 | + def truncate(self, point): |
| 212 | + if point.x >= self.maxx: |
| 213 | + point.x = self.maxx |
| 214 | + if point.y >= self.maxy: |
| 215 | + point.y = self.maxy |
| 216 | + if point.x < self.minx: |
| 217 | + point.x = self.minx |
| 218 | + if point.y < self.miny: |
| 219 | + point.y = self.miny |
| 220 | + |
| 221 | + def create_image_chaos(self, w, h, n): |
| 222 | + im = [[1] * h for i in range(w)] |
| 223 | + point = GVector((self.maxx + self.minx) / 2, |
| 224 | + (self.maxy + self.miny) / 2, 0) |
| 225 | + colored = 0 |
| 226 | + for _ in range(n): |
| 227 | + for i in range(5000): |
| 228 | + point = self.transform_point(point) |
| 229 | + x = (point.x - self.minx) / self.width * w |
| 230 | + y = (point.y - self.miny) / self.height * h |
| 231 | + x = int(x) |
| 232 | + y = int(y) |
| 233 | + if x == w: |
| 234 | + x -= 1 |
| 235 | + if y == h: |
| 236 | + y -= 1 |
| 237 | + im[x][h - y - 1] = 0 |
| 238 | + |
| 239 | + |
| 240 | +class Data: |
| 241 | + def __init__(self): |
| 242 | + self.c = None |
| 243 | + |
| 244 | +data = Data() |
| 245 | + |
| 246 | +def __setup__(iterations=10): |
| 247 | + splines = [ |
| 248 | + Spline([ |
| 249 | + GVector(1.597350, 3.304460, 0.000000), |
| 250 | + GVector(1.575810, 4.123260, 0.000000), |
| 251 | + GVector(1.313210, 5.288350, 0.000000), |
| 252 | + GVector(1.618900, 5.329910, 0.000000), |
| 253 | + GVector(2.889940, 5.502700, 0.000000), |
| 254 | + GVector(2.373060, 4.381830, 0.000000), |
| 255 | + GVector(1.662000, 4.360280, 0.000000)], |
| 256 | + 3, [0, 0, 0, 1, 1, 1, 2, 2, 2]), |
| 257 | + Spline([ |
| 258 | + GVector(2.804500, 4.017350, 0.000000), |
| 259 | + GVector(2.550500, 3.525230, 0.000000), |
| 260 | + GVector(1.979010, 2.620360, 0.000000), |
| 261 | + GVector(1.979010, 2.620360, 0.000000)], |
| 262 | + 3, [0, 0, 0, 1, 1, 1]), |
| 263 | + Spline([ |
| 264 | + GVector(2.001670, 4.011320, 0.000000), |
| 265 | + GVector(2.335040, 3.312830, 0.000000), |
| 266 | + GVector(2.366800, 3.233460, 0.000000), |
| 267 | + GVector(2.366800, 3.233460, 0.000000)], |
| 268 | + 3, [0, 0, 0, 1, 1, 1]) |
| 269 | + ] |
| 270 | + data.c = Chaosgame(splines, 0.25) |
| 271 | + |
| 272 | + |
| 273 | +def main(iterations=10, w=1000, h=1200): |
| 274 | + data.c.create_image_chaos(w, h, iterations) |
| 275 | + |
| 276 | +def __benchmark__(iterations=10): |
| 277 | + main(iterations) |
0 commit comments