topo.misc.patternfn

1 from __future__ import with_statement 2 """ 3 Family of two-dimensional functions indexed by x and y. 4 5 All functions are written to be valid both for scalar x and y, and for 6 numpy arrays of x and y (in which case the result is also an array); 7 the functions therefore have the same mathematical behaviour as numpy. 8 9 $Id: patternfn.py 11295 2010-07-27 14:21:39Z ceball $ 10 """ 11 __version__='$Revision: 11295 $' 12 13 14 15 from math import pi 16 17 from numpy.oldnumeric import where,maximum,cos,sqrt,divide,greater_equal,bitwise_xor,exp 18 from numpy.oldnumeric import arcsin,logical_and,logical_or,less,minimum 19 from numpy import seterr 20 21 from contextlib import contextmanager

22 23 @contextmanager 24 -def float_error_ignore():

25 """ 26 Many of the functions in this module use Gaussian smoothing, which 27 is based on a calculation like exp(divide(x*x,sigma)). When sigma 28 is zero the value of this expression should be zero at all points 29 in the plane, because such a Gaussian is infinitely small. 30 Obtaining the correct answer using finite-precision floating-point 31 array computations requires allowing infinite values to be 32 returned from divide(), and allowing exp() to underflow silently 33 to zero when given an infinite value. In numpy this is achieved 34 by using its seterr() function to disable divide-by-zero and 35 underflow warnings temporarily while these values are being 36 computed. 37 """ 38 oldsettings=seterr(divide='ignore',under='ignore') 39 yield 40 seterr(**oldsettings)

41

42 43 -def gaussian(x, y, xsigma, ysigma):

44 """ 45 Two-dimensional oriented Gaussian pattern (i.e., 2D version of a 46 bell curve, like a normal distribution but not necessarily summing 47 to 1.0). 48 """ 49 if xsigma==0.0 or ysigma==0.0: 50 return x*0.0 51 52 with float_error_ignore(): 53 x_w = divide(x,xsigma) 54 y_h = divide(y,ysigma) 55 return exp(-0.5*x_w*x_w + -0.5*y_h*y_h)

56

57 58 -def sigmoid(axis, slope):

59 """ 60 Sigmoid dividing axis into a positive and negative half, 61 with a smoothly sloping transition between them (controlled by the slope). 62 63 At default rotation, axis refers to the vertical (y) axis. 64 """ 65 with float_error_ignore(): 66 return (2.0 / (1.0 + exp(-2.0*slope*axis))) - 1.0 67

68 69 -def exponential(x, y, xscale, yscale):

70 """ 71 Two-dimensional oriented exponential decay pattern. 72 """ 73 if xscale==0.0 or yscale==0.0: 74 return x*0.0 75 76 with float_error_ignore(): 77 x_w = divide(x,xscale) 78 y_h = divide(y,yscale) 79 return exp(-sqrt(x_w*x_w+y_h*y_h))

80

81 82 -def gabor(x, y, xsigma, ysigma, frequency, phase):

83 """ 84 Gabor pattern (sine grating multiplied by a circular Gaussian). 85 """ 86 if xsigma==0.0 or ysigma==0.0: 87 return x*0.0 88 89 with float_error_ignore(): 90 x_w = divide(x,xsigma) 91 y_h = divide(y,ysigma) 92 p = exp(-0.5*x_w*x_w + -0.5*y_h*y_h) 93 return p * 0.5*cos(2*pi*frequency*y + phase)

94

95 96 # JABHACKALERT: Shouldn't this use 'size' instead of 'thickness', 97 # for consistency with the other patterns? Right now, it has a 98 # size parameter and ignores it, which is very confusing. I guess 99 # it's called thickness to match ring, but matching gaussian and disk 100 # is probably more important. 101 -def line(y, thickness, gaussian_width):

102 """ 103 Infinite-length line with a solid central region, then Gaussian fall-off at the edges. 104 """ 105 distance_from_line = abs(y) 106 gaussian_y_coord = distance_from_line - thickness/2.0 107 sigmasq = gaussian_width*gaussian_width 108 109 if sigmasq==0.0: 110 falloff = y*0.0 111 else: 112 with float_error_ignore(): 113 falloff = exp(divide(-gaussian_y_coord*gaussian_y_coord,2*sigmasq)) 114 115 return where(gaussian_y_coord<=0, 1.0, falloff)

116

117 118 -def disk(x, y, height, gaussian_width):

119 """ 120 Circular disk with Gaussian fall-off after the solid central region. 121 """ 122 disk_radius = height/2.0 123 124 distance_from_origin = sqrt(x**2+y**2) 125 distance_outside_disk = distance_from_origin - disk_radius 126 sigmasq = gaussian_width*gaussian_width 127 128 if sigmasq==0.0: 129 falloff = x*0.0 130 else: 131 with float_error_ignore(): 132 falloff = exp(divide(-distance_outside_disk*distance_outside_disk, 133 2*sigmasq)) 134 135 return where(distance_outside_disk<=0,1.0,falloff)

136

137 138 -def ring(x, y, height, thickness, gaussian_width):

139 """ 140 Circular ring (annulus) with Gaussian fall-off after the solid ring-shaped region. 141 """ 142 radius = height/2.0 143 half_thickness = thickness/2.0 144 145 distance_from_origin = sqrt(x**2+y**2) 146 distance_outside_outer_disk = distance_from_origin - radius - half_thickness 147 distance_inside_inner_disk = radius - half_thickness - distance_from_origin 148 149 ring = 1.0-bitwise_xor(greater_equal(distance_inside_inner_disk,0.0),greater_equal(distance_outside_outer_disk,0.0)) 150 151 sigmasq = gaussian_width*gaussian_width 152 153 if sigmasq==0.0: 154 inner_falloff = x*0.0 155 outer_falloff = x*0.0 156 else: 157 with float_error_ignore(): 158 inner_falloff = exp(divide(-distance_inside_inner_disk*distance_inside_inner_disk, 2.0*sigmasq)) 159 outer_falloff = exp(divide(-distance_outside_outer_disk*distance_outside_outer_disk, 2.0*sigmasq)) 160 161 return maximum(inner_falloff,maximum(outer_falloff,ring))

162

163 164 -def smooth_rectangle(x, y, rec_w, rec_h, gaussian_width_x, gaussian_width_y):

165 """ 166 Rectangle with a solid central region, then Gaussian fall-off at the edges. 167 """ 168 169 gaussian_x_coord = abs(x)-rec_w/2.0 170 gaussian_y_coord = abs(y)-rec_h/2.0 171 172 box_x=less(gaussian_x_coord,0.0) 173 box_y=less(gaussian_y_coord,0.0) 174 sigmasq_x=gaussian_width_x*gaussian_width_x 175 sigmasq_y=gaussian_width_y*gaussian_width_y 176 177 with float_error_ignore(): 178 falloff_x=x*0.0 if sigmasq_x==0.0 else \ 179 exp(divide(-gaussian_x_coord*gaussian_x_coord,2*sigmasq_x)) 180 falloff_y=y*0.0 if sigmasq_y==0.0 else \ 181 exp(divide(-gaussian_y_coord*gaussian_y_coord,2*sigmasq_y)) 182 183 return minimum(maximum(box_x,falloff_x), maximum(box_y,falloff_y))

184

185 186 187 -def arc_by_radian(x, y, height, radian_range, thickness, gaussian_width):

188 """ 189 Radial arc with Gaussian fall-off after the solid ring-shaped 190 region with the given thickness, with shape specified by the 191 (start,end) radian_range. 192 """ 193 194 # Create a circular ring (copied from the ring function) 195 radius = height/2.0 196 half_thickness = thickness/2.0 197 198 distance_from_origin = sqrt(x**2+y**2) 199 distance_outside_outer_disk = distance_from_origin - radius - half_thickness 200 distance_inside_inner_disk = radius - half_thickness - distance_from_origin 201 202 ring = 1.0-bitwise_xor(greater_equal(distance_inside_inner_disk,0.0),greater_equal(distance_outside_outer_disk,0.0)) 203 204 sigmasq = gaussian_width*gaussian_width 205 206 if sigmasq==0.0: 207 inner_falloff = x*0.0 208 outer_falloff = x*0.0 209 else: 210 with float_error_ignore(): 211 inner_falloff = exp(divide(-distance_inside_inner_disk*distance_inside_inner_disk, 2.0*sigmasq)) 212 outer_falloff = exp(divide(-distance_outside_outer_disk*distance_outside_outer_disk, 2.0*sigmasq)) 213 214 output_ring = maximum(inner_falloff,maximum(outer_falloff,ring)) 215 216 # Calculate radians (in 4 phases) and cut according to the set range) 217 218 # RZHACKALERT: 219 # Function float_error_ignore() cannot catch the exception when 220 # both dividend and divisor are 0.0, and when only divisor is 0.0 221 # it returns 'Inf' rather than 0.0. In x, y and 222 # distance_from_origin, only one point in distance_from_origin can 223 # be 0.0 (circle center) and in this point x and y must be 0.0 as 224 # well. So here is a hack to avoid the 'invalid value encountered 225 # in divide' error by turning 0.0 to 1e-5 in distance_from_origin. 226 distance_from_origin += where(distance_from_origin == 0.0, 1e-5, 0) 227 228 with float_error_ignore(): 229 sines = divide(y, distance_from_origin) 230 cosines = divide(x, distance_from_origin) 231 arcsines = arcsin(sines) 232 233 phase_1 = where(logical_and(sines >= 0, cosines >= 0), 2*pi-arcsines, 0) 234 phase_2 = where(logical_and(sines >= 0, cosines < 0), pi+arcsines, 0) 235 phase_3 = where(logical_and(sines < 0, cosines < 0), pi+arcsines, 0) 236 phase_4 = where(logical_and(sines < 0, cosines >= 0), -arcsines, 0) 237 arcsines = phase_1 + phase_2 + phase_3 + phase_4 238 239 if radian_range[0] <= radian_range[1]: 240 return where(logical_and(arcsines >= radian_range[0], arcsines <= radian_range[1]), 241 output_ring, 0.0) 242 else: 243 return where(logical_or(arcsines >= radian_range[0], arcsines <= radian_range[1]), 244 output_ring, 0.0)

245

246 247 -def arc_by_center(x, y, arc_box, constant_length, thickness, gaussian_width):

248 """ 249 Arc with Gaussian fall-off after the solid ring-shaped region and specified 250 by point of tangency (x and y) and arc width and height. 251 252 This function calculates the start and end radian from the given width and 253 height, and then calls arc_by_radian function to draw the curve. 254 """ 255 256 arc_w=arc_box[0] 257 arc_h=abs(arc_box[1]) 258 259 if arc_w==0.0: # arc_w=0, don't draw anything 260 radius=0.0 261 angles=(0.0,0.0) 262 elif arc_h==0.0: # draw a horizontal line, width=arc_w 263 return smooth_rectangle(x, y, arc_w, thickness, 0.0, gaussian_width) 264 else: 265 if constant_length: 266 curvature=arc_h/arc_w 267 radius=arc_w/(2*pi*curvature) 268 angle=curvature*(2*pi)/2.0 269 else: # constant width 270 radius=arc_h/2.0+arc_w**2.0/(8*arc_h) 271 angle=arcsin(arc_w/2.0/radius) 272 if arc_box[1]<0: # convex shape 273 y=y+radius 274 angles=(3.0/2.0*pi-angle, 3.0/2.0*pi+angle) 275 else: # concave shape 276 y=y-radius 277 angles=(pi/2.0-angle, pi/2.0+angle) 278 279 return arc_by_radian(x, y, radius*2.0, angles, thickness, gaussian_width)

280

Source Code for Module topo.misc.patternfn