|
| 1 | + |
| 2 | +# coding: utf-8 |
| 3 | + |
| 4 | +# In[1]: |
| 5 | + |
| 6 | + |
| 7 | +#another fractal script via recursion |
| 8 | +#this script heavily involves analytic geometry |
| 9 | +#a good material to help you rewind high school geometry |
| 10 | +# https://youthconway.com/handouts/analyticgeo.pdf |
| 11 | +import matplotlib.pyplot as plt |
| 12 | + |
| 13 | + |
| 14 | +# In[2]: |
| 15 | + |
| 16 | + |
| 17 | +#simple solution to get coefficients of the equation |
| 18 | +def get_line_params(x1,y1,x2,y2): |
| 19 | + |
| 20 | + slope=(y1-y2)/(x1-x2) |
| 21 | + intercept=y1-slope*x1 |
| 22 | + |
| 23 | + return slope,intercept |
| 24 | + |
| 25 | + |
| 26 | +# In[3]: |
| 27 | + |
| 28 | + |
| 29 | +#compute euclidean distance |
| 30 | +def euclidean_distance(point1,point2): |
| 31 | + return ((point1[0]-point2[0])**2+(point1[1]-point2[1])**2)**0.5 |
| 32 | + |
| 33 | + |
| 34 | +# In[4]: |
| 35 | + |
| 36 | + |
| 37 | +#standard solution to quadratic equation |
| 38 | +def solve_quadratic_equation(A,B,C): |
| 39 | + x1=(-B+(B**2-4*A*C)**0.5)/(2*A) |
| 40 | + x2=(-B-(B**2-4*A*C)**0.5)/(2*A) |
| 41 | + return [x1,x2] |
| 42 | + |
| 43 | + |
| 44 | +# In[5]: |
| 45 | + |
| 46 | + |
| 47 | +#analytic geometry to compute target datapoints |
| 48 | +def get_datapoint(pivot,measure,length,direction='inner'): |
| 49 | + |
| 50 | + #for undefined slope |
| 51 | + if pivot[0]==measure[0]: |
| 52 | + y1=pivot[1]+length |
| 53 | + y2=pivot[1]-length |
| 54 | + x1=pivot[0] |
| 55 | + x2=pivot[0] |
| 56 | + |
| 57 | + #for general cases |
| 58 | + else: |
| 59 | + |
| 60 | + #get line equation |
| 61 | + slope,intercept=get_line_params(pivot[0],pivot[1], |
| 62 | + measure[0],measure[1],) |
| 63 | + |
| 64 | + #solve quadratic equation |
| 65 | + A=1 |
| 66 | + B=-2*pivot[0] |
| 67 | + C=pivot[0]**2-length**2/(slope**2+1) |
| 68 | + x1,x2=solve_quadratic_equation(A,B,C) |
| 69 | + |
| 70 | + #get y from line equation |
| 71 | + y1=slope*x1+intercept |
| 72 | + y2=slope*x2+intercept |
| 73 | + |
| 74 | + if direction=='inner': |
| 75 | + |
| 76 | + #take the one between pivot and measure points |
| 77 | + datapoint=min([(x1,y1),(x2,y2)], |
| 78 | + key=lambda x:euclidean_distance(x,measure)) |
| 79 | + else: |
| 80 | + |
| 81 | + #take the one farther away from measure points |
| 82 | + datapoint=max([(x1,y1),(x2,y2)], |
| 83 | + key=lambda x:euclidean_distance(x,measure)) |
| 84 | + |
| 85 | + return datapoint |
| 86 | + |
| 87 | + |
| 88 | +# In[6]: |
| 89 | + |
| 90 | + |
| 91 | +#recursively compute the coordinates of koch curve data points |
| 92 | +#to effectively connect the data points,the best choice is to use turtle |
| 93 | +#it would be too difficult to connect the dots via analytic geometry |
| 94 | +def koch_snowflake(base1,base2,base3,n): |
| 95 | + |
| 96 | + #base case |
| 97 | + if n==0: |
| 98 | + return |
| 99 | + |
| 100 | + else: |
| 101 | + |
| 102 | + #find mid point |
| 103 | + #midpoint between base1 and base2 has to satisfy two conditions |
| 104 | + #it has to be on the same line as base1 and base2 |
| 105 | + #assume this line follows y=kx+b |
| 106 | + #the midpoint is (x,kx+b) |
| 107 | + #base1 is (α,kα+b),base2 is (δ,kδ+b) |
| 108 | + #the euclidean distance between midpoint and base1 should be |
| 109 | + #half of the euclidean distance between base1 and base2 |
| 110 | + #(x-α)**2+(kx+b-kα-b)**2=((α-δ)**2+(kα+b-kδ-b)**2)/4 |
| 111 | + #apart from x,everything else in the equation is constant |
| 112 | + #this forms a simple quadratic solution to get two roots |
| 113 | + #one root would be between base1 and base2 which yields midpoint |
| 114 | + #and the other would be farther away from base2 |
| 115 | + #this function solves the equation via (-B+(B**2-4*A*C)**0.5)/(2*A) |
| 116 | + #alternatively,you can use scipy.optimize.root |
| 117 | + #the caveat is it does not offer both roots |
| 118 | + #a wrong initial guess could take you to the wrong root |
| 119 | + midpoint=get_datapoint(base1,base2,euclidean_distance(base1,base2)/2) |
| 120 | + |
| 121 | + #compute the top point of a triangle |
| 122 | + #the computation is similar to midpoint |
| 123 | + #the euclidean distance between triangle_top and midpoint should be |
| 124 | + #one third of the distance between base3 and midpoint |
| 125 | + triangle_top=get_datapoint(midpoint,base3, |
| 126 | + euclidean_distance(midpoint,base3)/3, |
| 127 | + direction='outer') |
| 128 | + |
| 129 | + #two segment points divide a line into three equal parts |
| 130 | + #the computation is almost the same as midpoint |
| 131 | + #the euclidean distance between segment1 and base1 |
| 132 | + #should be one third of the euclidean distance between base2 and base1 |
| 133 | + segment1=get_datapoint(base1,base2,euclidean_distance(base1,base2)/3) |
| 134 | + segment2=get_datapoint(base2,base1,euclidean_distance(base1,base2)/3) |
| 135 | + |
| 136 | + #compute the nearest segment point of the neighboring line |
| 137 | + segment_side_1=get_datapoint(base1,base3,euclidean_distance(base1,base3)/3) |
| 138 | + segment_side_2=get_datapoint(base2,base3,euclidean_distance(base2,base3)/3) |
| 139 | + |
| 140 | + #recursion |
| 141 | + yield [segment1,segment2,triangle_top] |
| 142 | + yield from koch_snowflake(base1,segment1,segment_side_1,n-1) |
| 143 | + yield from koch_snowflake(segment1,triangle_top,segment2,n-1) |
| 144 | + yield from koch_snowflake(triangle_top,segment2,segment1,n-1) |
| 145 | + yield from koch_snowflake(segment2,base2,segment_side_2,n-1) |
| 146 | + |
| 147 | + |
| 148 | +# In[7]: |
| 149 | + |
| 150 | + |
| 151 | +#set data points |
| 152 | +point1=(0,0) |
| 153 | +point2=(3,0) |
| 154 | +point3=(3/2,3/2*(3**0.5)) |
| 155 | + |
| 156 | +#due to python floating point arithmetic |
| 157 | +#a lot of data points could go wrong during the calculation |
| 158 | +#unfortunately there is no panacea to this malaise |
| 159 | +#unless we plan to use decimal package all the time |
| 160 | +#when the depth of snowflake reaches 1 |
| 161 | +#one of the data points reach -1.1102230246251565e-16 on x axis |
| 162 | +#when the depth of snowflake reaches 6 |
| 163 | +#we end up with complex root and things go wrong |
| 164 | +n=4 |
| 165 | + |
| 166 | + |
| 167 | +# In[8]: |
| 168 | + |
| 169 | + |
| 170 | +#collect coordinates |
| 171 | +arr=list(koch_snowflake(point1,point2,point3,n))+list( |
| 172 | + koch_snowflake(point1,point3,point2,n))+list( |
| 173 | + koch_snowflake(point2,point3,point1,n))+[(point1,point2,point3)] |
| 174 | +coordinates=[j for i in arr for j in i] |
| 175 | + |
| 176 | + |
| 177 | +# In[9]: |
| 178 | + |
| 179 | + |
| 180 | +#viz |
| 181 | +#visually u can tell some of the data points are miscalculated |
| 182 | +#purely caused by floating point arithmetic |
| 183 | +ax=plt.figure(figsize=(5,5)) |
| 184 | +plt.scatter([i[0] for i in coordinates], |
| 185 | + [i[1] for i in coordinates],s=0.5) |
| 186 | +plt.axis('off') |
| 187 | +plt.show() |
0 commit comments