draw3d.py 3.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107
  1. import numpy as np
  2. from scipy.special import elliprg
  3. def ellipsoid(a, b, c, spacing=(1.0, 1.0, 1.0), levelset=False):
  4. """Generate ellipsoid for given semi-axis lengths.
  5. The respective semi-axis lengths are given along three dimensions in
  6. Cartesian coordinates. Each dimension may use a different grid spacing.
  7. Parameters
  8. ----------
  9. a : float
  10. Length of semi-axis along x-axis.
  11. b : float
  12. Length of semi-axis along y-axis.
  13. c : float
  14. Length of semi-axis along z-axis.
  15. spacing : 3-tuple of floats
  16. Grid spacing in three spatial dimensions.
  17. levelset : bool
  18. If True, returns the level set for this ellipsoid (signed level
  19. set about zero, with positive denoting interior) as np.float64.
  20. False returns a binarized version of said level set.
  21. Returns
  22. -------
  23. ellipsoid : (M, N, P) array
  24. Ellipsoid centered in a correctly sized array for given `spacing`.
  25. Boolean dtype unless `levelset=True`, in which case a float array is
  26. returned with the level set above 0.0 representing the ellipsoid.
  27. """
  28. if (a <= 0) or (b <= 0) or (c <= 0):
  29. raise ValueError('Parameters a, b, and c must all be > 0')
  30. offset = np.r_[1, 1, 1] * np.r_[spacing]
  31. # Calculate limits, and ensure output volume is odd & symmetric
  32. low = np.ceil(-np.r_[a, b, c] - offset)
  33. high = np.floor(np.r_[a, b, c] + offset + 1)
  34. for dim in range(3):
  35. if (high[dim] - low[dim]) % 2 == 0:
  36. low[dim] -= 1
  37. num = np.arange(low[dim], high[dim], spacing[dim])
  38. if 0 not in num:
  39. low[dim] -= np.max(num[num < 0])
  40. # Generate (anisotropic) spatial grid
  41. x, y, z = np.mgrid[
  42. low[0] : high[0] : spacing[0],
  43. low[1] : high[1] : spacing[1],
  44. low[2] : high[2] : spacing[2],
  45. ]
  46. if not levelset:
  47. arr = ((x / float(a)) ** 2 + (y / float(b)) ** 2 + (z / float(c)) ** 2) <= 1
  48. else:
  49. arr = ((x / float(a)) ** 2 + (y / float(b)) ** 2 + (z / float(c)) ** 2) - 1
  50. return arr
  51. def ellipsoid_stats(a, b, c):
  52. """Calculate analytical volume and surface area of an ellipsoid.
  53. The surface area of an ellipsoid is given by
  54. .. math:: S=4\\pi b c R_G\\!\\left(1, \\frac{a^2}{b^2}, \\frac{a^2}{c^2}\\right)
  55. where :math:`R_G` is Carlson's completely symmetric elliptic integral of
  56. the second kind [1]_. The latter is implemented as
  57. :py:func:`scipy.special.elliprg`.
  58. Parameters
  59. ----------
  60. a : float
  61. Length of semi-axis along x-axis.
  62. b : float
  63. Length of semi-axis along y-axis.
  64. c : float
  65. Length of semi-axis along z-axis.
  66. Returns
  67. -------
  68. vol : float
  69. Calculated volume of ellipsoid.
  70. surf : float
  71. Calculated surface area of ellipsoid.
  72. References
  73. ----------
  74. .. [1] Paul Masson (2020). Surface Area of an Ellipsoid.
  75. https://analyticphysics.com/Mathematical%20Methods/Surface%20Area%20of%20an%20Ellipsoid.htm
  76. """
  77. if (a <= 0) or (b <= 0) or (c <= 0):
  78. raise ValueError('Parameters a, b, and c must all be > 0')
  79. # Volume
  80. vol = 4 / 3.0 * np.pi * a * b * c
  81. # Surface area
  82. surf = 3 * vol * elliprg(1 / a**2, 1 / b**2, 1 / c**2)
  83. return vol, surf