Use hypot method to avoid over/underflow errors

Use hypot instead of manual calculation
This commit is contained in:
Ogi Moore 2021-04-17 22:38:21 -07:00
parent b0a3849960
commit f8cefa6284
10 changed files with 27 additions and 40 deletions

View File

@ -28,7 +28,7 @@ def psi(i, j, k, offset=(25, 25, 50)):
x = i-offset[0] x = i-offset[0]
y = j-offset[1] y = j-offset[1]
z = k-offset[2] z = k-offset[2]
th = np.arctan2(z, np.sqrt(x**2+y**2)) th = np.arctan2(z, np.hypot(x, y))
phi = np.arctan2(y, x) phi = np.arctan2(y, x)
r = np.sqrt(x**2 + y**2 + z **2) r = np.sqrt(x**2 + y**2 + z **2)
a0 = 1 a0 = 1

View File

@ -30,14 +30,14 @@ gz.translate(0, 0, -10)
w.addItem(gz) w.addItem(gz)
def fn(x, y): def fn(x, y):
return np.cos((x**2 + y**2)**0.5) return np.cos(np.hypot(x, y))
n = 51 n = 51
y = np.linspace(-10,10,n) y = np.linspace(-10,10,n)
x = np.linspace(-10,10,100) x = np.linspace(-10,10,100)
for i in range(n): for i in range(n):
yi = np.array([y[i]]*100) yi = np.array([y[i]]*100)
d = (x**2 + yi**2)**0.5 d = np.hypot(x, yi)
z = 10 * np.cos(d) / (d+1) z = 10 * np.cos(d) / (d+1)
pts = np.vstack([x,yi,z]).transpose() pts = np.vstack([x,yi,z]).transpose()
plt = gl.GLLinePlotItem(pos=pts, color=pg.glColor((i,n*1.3)), width=(i+1)/10., antialias=True) plt = gl.GLLinePlotItem(pos=pts, color=pg.glColor((i,n*1.3)), width=(i+1)/10., antialias=True)

View File

@ -29,17 +29,14 @@ def psi(i, j, k, offset=(50,50,100)):
x = i-offset[0] x = i-offset[0]
y = j-offset[1] y = j-offset[1]
z = k-offset[2] z = k-offset[2]
th = np.arctan2(z, (x**2+y**2)**0.5) th = np.arctan2(z, np.hypot(x, y))
phi = np.arctan2(y, x) phi = np.arctan2(y, x)
r = (x**2 + y**2 + z **2)**0.5 r = np.sqrt(x**2 + y**2 + z **2)
a0 = 2 a0 = 2
#ps = (1./81.) * (2./np.pi)**0.5 * (1./a0)**(3/2) * (6 - r/a0) * (r/a0) * np.exp(-r/(3*a0)) * np.cos(th) #ps = (1./81.) * (2./np.pi)**0.5 * (1./a0)**(3/2) * (6 - r/a0) * (r/a0) * np.exp(-r/(3*a0)) * np.cos(th)
ps = (1./81.) * 1./(6.*np.pi)**0.5 * (1./a0)**(3/2) * (r/a0)**2 * np.exp(-r/(3*a0)) * (3 * np.cos(th)**2 - 1) ps = (1./81.) * 1./(6.*np.pi)**0.5 * (1./a0)**(3/2) * (r/a0)**2 * np.exp(-r/(3*a0)) * (3 * np.cos(th)**2 - 1)
return ps return ps
#return ((1./81.) * (1./np.pi)**0.5 * (1./a0)**(3/2) * (r/a0)**2 * (r/a0) * np.exp(-r/(3*a0)) * np.sin(th) * np.cos(th) * np.exp(2 * 1j * phi))**2
data = np.fromfunction(psi, (100,100,200)) data = np.fromfunction(psi, (100,100,200))
positive = np.log(np.clip(data, 0, data.max())**2) positive = np.log(np.clip(data, 0, data.max())**2)

View File

@ -1,5 +1,5 @@
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
from math import atan2, asin, sin, cos, sqrt, pi from math import atan2, asin, sin, cos, sqrt, pi, hypot
import pyqtgraph as pg import pyqtgraph as pg
from pyqtgraph.Qt import QtGui, QtCore from pyqtgraph.Qt import QtGui, QtCore
import numpy as np import numpy as np
@ -417,7 +417,7 @@ class CircleSurface(pg.GraphicsObject):
## find intersection of circle and line (quadratic formula) ## find intersection of circle and line (quadratic formula)
dx = dir[0] dx = dir[0]
dy = dir[1] dy = dir[1]
dr = (dx**2 + dy**2) ** 0.5 dr = hypot(dx, dy) # length
D = p[0] * (p[1]+dy) - (p[0]+dx) * p[1] D = p[0] * (p[1]+dy) - (p[0]+dx) * p[1]
idr2 = 1.0 / dr**2 idr2 = 1.0 / dr**2
disc = r**2 * dr**2 - D**2 disc = r**2 * dr**2 - D**2
@ -429,7 +429,6 @@ class CircleSurface(pg.GraphicsObject):
else: else:
sgn = 1 sgn = 1
br = self.path.boundingRect() br = self.path.boundingRect()
x1 = (D*dy + sgn*dx*disc2) * idr2 x1 = (D*dy + sgn*dx*disc2) * idr2
y1 = (-D*dx + abs(dy)*disc2) * idr2 y1 = (-D*dx + abs(dy)*disc2) * idr2

View File

@ -6,7 +6,7 @@ Distributed under MIT/X11 license. See license.txt for more information.
""" """
from .Qt import QtCore from .Qt import QtCore
from math import sin, acos, atan2, inf, pi from math import sin, acos, atan2, inf, pi, hypot
def clip(x, mn, mx): def clip(x, mn, mx):
if x > mx: if x > mx:
@ -93,25 +93,13 @@ class Point(QtCore.QPointF):
return self._math_('__pow__', a) return self._math_('__pow__', a)
def _math_(self, op, x): def _math_(self, op, x):
#print "point math:", op
#try:
#fn = getattr(QtCore.QPointF, op)
#pt = fn(self, x)
#print fn, pt, self, x
#return Point(pt)
#except AttributeError:
x = Point(x) x = Point(x)
return Point(getattr(self[0], op)(x[0]), getattr(self[1], op)(x[1])) return Point(getattr(self[0], op)(x[0]), getattr(self[1], op)(x[1]))
def length(self): def length(self):
"""Returns the vector length of this Point.""" """Returns the vector length of this Point."""
try: return hypot(self[0], self[1]) # length
return (self[0]**2 + self[1]**2) ** 0.5
except OverflowError:
try:
return self[1] / sin(atan2(self[1], self[0]))
except OverflowError:
return inf
def norm(self): def norm(self):
"""Returns a vector in the same direction with unit length.""" """Returns a vector in the same direction with unit length."""

View File

@ -1,3 +1,4 @@
from math import hypot
from ..Qt import QtGui, QtCore from ..Qt import QtGui, QtCore
from .. import functions as fn from .. import functions as fn
import numpy as np import numpy as np
@ -135,7 +136,7 @@ class ArrowItem(QtGui.QGraphicsPathItem):
pad = 0 pad = 0
if self.opts['pxMode']: if self.opts['pxMode']:
br = self.boundingRect() br = self.boundingRect()
pad += (br.width()**2 + br.height()**2) ** 0.5 pad += hypot(br.width(), br.height())
pen = self.pen() pen = self.pen()
if pen.isCosmetic(): if pen.isCosmetic():
pad += max(1, pen.width()) * 0.7072 pad += max(1, pen.width()) * 0.7072

View File

@ -1,4 +1,5 @@
import warnings import warnings
from math import hypot
from collections import OrderedDict from collections import OrderedDict
from functools import reduce from functools import reduce
from ..Qt import QtGui, QtCore, isQObjectAlive from ..Qt import QtGui, QtCore, isQObjectAlive
@ -308,7 +309,7 @@ class GraphicsItem(object):
v = self.pixelVectors() v = self.pixelVectors()
if v == (None, None): if v == (None, None):
return None, None return None, None
return (v[0].x()**2+v[0].y()**2)**0.5, (v[1].x()**2+v[1].y()**2)**0.5 return (hypot(v[0].x(), v[0].y()), hypot(v[1].x(), v[1].y())) # lengths
def pixelWidth(self): def pixelWidth(self):
## deprecated ## deprecated

View File

@ -17,7 +17,7 @@ import numpy as np
#from numpy.linalg import norm #from numpy.linalg import norm
from ..Point import * from ..Point import *
from ..SRTTransform import SRTTransform from ..SRTTransform import SRTTransform
from math import atan2, cos, sin, pi, sqrt from math import atan2, cos, sin, pi, sqrt, hypot
from .. import functions as fn from .. import functions as fn
from .GraphicsObject import GraphicsObject from .GraphicsObject import GraphicsObject
from .UIGraphicsItem import UIGraphicsItem from .UIGraphicsItem import UIGraphicsItem
@ -1247,8 +1247,8 @@ class ROI(GraphicsObject):
vx = img.mapToData(self.mapToItem(img, QtCore.QPointF(1, 0))) - origin vx = img.mapToData(self.mapToItem(img, QtCore.QPointF(1, 0))) - origin
vy = img.mapToData(self.mapToItem(img, QtCore.QPointF(0, 1))) - origin vy = img.mapToData(self.mapToItem(img, QtCore.QPointF(0, 1))) - origin
lvx = sqrt(vx.x()**2 + vx.y()**2) lvx = hypot(vx.x(), vx.y()) # length
lvy = sqrt(vy.x()**2 + vy.y()**2) lvy = hypot(vy.x(), vy.y()) # length
##img.width is number of pixels, not width of item. ##img.width is number of pixels, not width of item.
##need pxWidth and pxHeight instead of pxLen ? ##need pxWidth and pxHeight instead of pxLen ?
sx = 1.0 / lvx sx = 1.0 / lvx
@ -1887,7 +1887,7 @@ class EllipseROI(ROI):
h = arr.shape[axes[1]] h = arr.shape[axes[1]]
## generate an ellipsoidal mask ## generate an ellipsoidal mask
mask = np.fromfunction(lambda x,y: (((x+0.5)/(w/2.)-1)**2+ ((y+0.5)/(h/2.)-1)**2)**0.5 < 1, (w, h)) mask = np.fromfunction(lambda x,y: np.hypot(((x+0.5)/(w/2.)-1), ((y+0.5)/(h/2.)-1)) < 1, (w, h))
# reshape to match array axes # reshape to match array axes
if axes[0] > axes[1]: if axes[0] > axes[1]:

View File

@ -3,7 +3,8 @@ from OpenGL.GL import *
from OpenGL.arrays import vbo from OpenGL.arrays import vbo
from .. GLGraphicsItem import GLGraphicsItem from .. GLGraphicsItem import GLGraphicsItem
from .. import shaders from .. import shaders
from ... import QtGui from ...functions import clip_array
from ...Qt import QtGui
import numpy as np import numpy as np
__all__ = ['GLScatterPlotItem'] __all__ = ['GLScatterPlotItem']
@ -62,12 +63,11 @@ class GLScatterPlotItem(GLGraphicsItem):
## Generate texture for rendering points ## Generate texture for rendering points
w = 64 w = 64
def fn(x,y): def fn(x,y):
r = ((x-(w-1)/2.)**2 + (y-(w-1)/2.)**2) ** 0.5 r = np.hypot((x-(w-1)/2.), (y-(w-1)/2.))
return 255 * (w/2. - np.clip(r, w/2.-1.0, w/2.)) return 255 * (w/2. - clip_array(r, w/2.-1.0, w/2.))
pData = np.empty((w, w, 4)) pData = np.empty((w, w, 4))
pData[:] = 255 pData[:] = 255
pData[:,:,3] = np.fromfunction(fn, pData.shape[:2]) pData[:,:,3] = np.fromfunction(fn, pData.shape[:2])
#print pData.shape, pData.min(), pData.max()
pData = pData.astype(np.ubyte) pData = pData.astype(np.ubyte)
if getattr(self, "pointTexture", None) is None: if getattr(self, "pointTexture", None) is None:

View File

@ -1,4 +1,5 @@
from ..Qt import QtGui, QtCore from math import hypot
from ..Qt import QtGui, QtCore, mkQApp
__all__ = ['JoystickButton'] __all__ = ['JoystickButton']
@ -41,7 +42,7 @@ class JoystickButton(QtGui.QPushButton):
def setState(self, *xy): def setState(self, *xy):
xy = list(xy) xy = list(xy)
d = (xy[0]**2 + xy[1]**2)**0.5 d = hypot(xy[0], xy[1]) # length
nxy = [0, 0] nxy = [0, 0]
for i in [0,1]: for i in [0,1]:
if xy[i] == 0: if xy[i] == 0:
@ -84,7 +85,7 @@ class JoystickButton(QtGui.QPushButton):
if __name__ == '__main__': if __name__ == '__main__':
app = pg.mkQApp() app = mkQApp()
w = QtGui.QMainWindow() w = QtGui.QMainWindow()
b = JoystickButton() b = JoystickButton()
w.setCentralWidget(b) w.setCentralWidget(b)