我一直在尝试用Python编写自己的物理引擎，作为物理和编程方面的练习。我从这里的教程开始。一切进展顺利，但后来我发现了thomas jakobsen的文章《Advanced character physics》，其中介绍了如何使用Verlet积分进行模拟。
我一直在尝试使用verlet集成编写我自己的基本物理模拟器，但事实证明比我最初预期的要稍微困难一些。我在外面浏览一些程序的例子，偶然发现了这个用Python写的，我还发现了这个使用Processing的教程。
Processing版本给我印象最深的是它的运行速度。单是布料就有2400个不同的点被模拟，这还不包括身体。
python的例子只使用了256个粒子作为布料，并且它以每秒30帧的速度运行。我尝试将粒子的数量增加到2401(它必须是方形的程序才能工作)，它以大约3 fps的速度运行。
这两种方法都是将粒子对象的实例存储在列表中，然后遍历列表，调用每个粒子的"更新位置"方法。作为一个例子，这是Processing草图中计算每个粒子的新位置的部分代码:
for (int i = 0; i < pointmasses.size(); i++) {
PointMass pointmass = (PointMass) pointmasses.get(i);
pointmass.updateInteractions();
pointmass.updatePhysics(fixedDeltaTimeSeconds);
}
EDIT:这是我之前链接的python版本的代码:
"""
verletCloth01.py
Eric Pavey - 2010-07-03 - www.akeric.com
Riding on the shoulders of giants.
I wanted to learn now to do 'verlet cloth' in PythonPygame. I first ran across
this post source:
http://forums.overclockers.com.au/showthread.php?t=870396
http://dl.dropbox.com/u/3240460/cloth5.py
Which pointed to some good reference, that was a dead link. After some searching,
I found it here:
http://www.gpgstudy.com/gpgiki/GDC%202001%3A%20Advanced%20Character%20Physics
Which is a 2001 SIGGRAPH paper by Thomas Jakobsen called:
"GDC 2001: Advanced Characer Physics".
This code is a PythonPygame interpretation of that 2001 Siggraph paper. I did
borrow some code from 'domlebo's source code, it was a great starting point. But
I'd like to think I put my own flavor on it.
"""
#--------------
# Imports & Initis
import sys
from math import sqrt
# Vec2D comes from here: http://pygame.org/wiki/2DVectorClass
from vec2d import Vec2d
import pygame
from pygame.locals import *
pygame.init()
#--------------
# Constants
TITLE = "verletCloth01"
WIDTH = 600
HEIGHT = 600
FRAMERATE = 60
# How many iterations to run on our constraints per frame?
# This will 'tighten' the cloth, but slow the sim.
ITERATE = 2
GRAVITY = Vec2d(0.0,0.05)
TSTEP = 2.8
# How many pixels to position between each particle?
PSTEP = int(WIDTH*.03)
# Offset in pixels from the top left of screen to position grid:
OFFSET = int(.25*WIDTH)
#-------------
# Define helper functions, classes
class Particle(object):
"""
Stores position, previous position, and where it is in the grid.
"""
def __init__(self, screen, currentPos, gridIndex):
# Current Position : m_x
self.currentPos = Vec2d(currentPos)
# Index [x][y] of Where it lives in the grid
self.gridIndex = gridIndex
# Previous Position : m_oldx
self.oldPos = Vec2d(currentPos)
# Force accumulators : m_a
self.forces = GRAVITY
# Should the particle be locked at its current position?
self.locked = False
self.followMouse = False
self.colorUnlocked = Color('white')
self.colorLocked = Color('green')
self.screen = screen
def __str__(self):
return "Particle <%s, %s>"%(self.gridIndex[0], self.gridIndex[1])
def draw(self):
# Draw a circle at the given Particle.
screenPos = (self.currentPos[0], self.currentPos[1])
if self.locked:
pygame.draw.circle(self.screen, self.colorLocked, (int(screenPos[0]),
int(screenPos[1])), 4, 0)
else:
pygame.draw.circle(self.screen, self.colorUnlocked, (int(screenPos[0]),
int(screenPos[1])), 1, 0)
class Constraint(object):
"""
Stores 'constraint' data between two Particle objects. Stores this data
before the sim runs, to speed sim and draw operations.
"""
def __init__(self, screen, particles):
self.particles = sorted(particles)
# Calculate restlength as the initial distance between the two particles:
self.restLength = sqrt(abs(pow(self.particles[1].currentPos.x -
self.particles[0].currentPos.x, 2) +
pow(self.particles[1].currentPos.y -
self.particles[0].currentPos.y, 2)))
self.screen = screen
self.color = Color('red')
def __str__(self):
return "Constraint <%s, %s>"%(self.particles[0], self.particles[1])
def draw(self):
# Draw line between the two particles.
p1 = self.particles[0]
p2 = self.particles[1]
p1pos = (p1.currentPos[0],
p1.currentPos[1])
p2pos = (p2.currentPos[0],
p2.currentPos[1])
pygame.draw.aaline(self.screen, self.color,
(p1pos[0], p1pos[1]), (p2pos[0], p2pos[1]), 1)
class Grid(object):
"""
Stores a grid of Particle objects. Emulates a 2d container object. Particle
objects can be indexed by position:
grid = Grid()
particle = g[2][4]
"""
def __init__(self, screen, rows, columns, step, offset):
self.screen = screen
self.rows = rows
self.columns = columns
self.step = step
self.offset = offset
# Make our internal grid:
# _grid is a list of sublists.
# Each sublist is a 'column'.
# Each column holds a particle object per row:
# _grid =
# [[p00, [p10, [etc,
# p01, p11,
# etc], etc], ]]
self._grid = []
for x in range(columns):
self._grid.append([])
for y in range(rows):
currentPos = (x*self.step+self.offset, y*self.step+self.offset)
self._grid[x].append(Particle(self.screen, currentPos, (x,y)))
def getNeighbors(self, gridIndex):
"""
return a list of all neighbor particles to the particle at the given gridIndex:
gridIndex = [x,x] : The particle index we're polling
"""
possNeighbors = []
possNeighbors.append([gridIndex[0]-1, gridIndex[1]])
possNeighbors.append([gridIndex[0], gridIndex[1]-1])
possNeighbors.append([gridIndex[0]+1, gridIndex[1]])
possNeighbors.append([gridIndex[0], gridIndex[1]+1])
neigh = []
for coord in possNeighbors:
if (coord[0] < 0) | (coord[0] > self.rows-1):
pass
elif (coord[1] < 0) | (coord[1] > self.columns-1):
pass
else:
neigh.append(coord)
finalNeighbors = []
for point in neigh:
finalNeighbors.append((point[0], point[1]))
return finalNeighbors
#--------------------------
# Implement Container Type:
def __len__(self):
return len(self.rows * self.columns)
def __getitem__(self, key):
return self._grid[key]
def __setitem__(self, key, value):
self._grid[key] = value
#def __delitem__(self, key):
#del(self._grid[key])
def __iter__(self):
for x in self._grid:
for y in x:
yield y
def __contains__(self, item):
for x in self._grid:
for y in x:
if y is item:
return True
return False
class ParticleSystem(Grid):
"""
Implements the verlet particles physics on the encapsulated Grid object.
"""
def __init__(self, screen, rows=49, columns=49, step=PSTEP, offset=OFFSET):
super(ParticleSystem, self).__init__(screen, rows, columns, step, offset)
# Generate our list of Constraint objects. One is generated between
# every particle connection.
self.constraints = []
for p in self:
neighborIndices = self.getNeighbors(p.gridIndex)
for ni in neighborIndices:
# Get the neighbor Particle from the index:
n = self[ni[0]][ni[1]]
# Let's not add duplicate Constraints, which would be easy to do!
new = True
for con in self.constraints:
if n in con.particles and p in con.particles:
new = False
if new:
self.constraints.append( Constraint(self.screen, (p,n)) )
# Lock our top left and right particles by default:
self[0][0].locked = True
self[1][0].locked = True
self[-2][0].locked = True
self[-1][0].locked = True
def verlet(self):
# Verlet integration step:
for p in self:
if not p.locked:
# make a copy of our current position
temp = Vec2d(p.currentPos)
p.currentPos += p.currentPos - p.oldPos + p.forces * TSTEP**2
p.oldPos = temp
elif p.followMouse:
temp = Vec2d(p.currentPos)
p.currentPos = Vec2d(pygame.mouse.get_pos())
p.oldPos = temp
def satisfyConstraints(self):
# Keep particles together:
for c in self.constraints:
delta = c.particles[0].currentPos - c.particles[1].currentPos
deltaLength = sqrt(delta.dot(delta))
try:
# You can get a ZeroDivisionError here once, so let's catch it.
# I think it's when particles sit on top of one another due to
# being locked.
diff = (deltaLength-c.restLength)/deltaLength
if not c.particles[0].locked:
c.particles[0].currentPos -= delta*0.5*diff
if not c.particles[1].locked:
c.particles[1].currentPos += delta*0.5*diff
except ZeroDivisionError:
pass
def accumulateForces(self):
# This doesn't do much right now, other than constantly reset the
# particles 'forces' to be 'gravity'. But this is where you'd implement
# other things, like drag, wind, etc.
for p in self:
p.forces = GRAVITY
def timeStep(self):
# This executes the whole shebang:
self.accumulateForces()
self.verlet()
for i in range(ITERATE):
self.satisfyConstraints()
def draw(self):
"""
Draw constraint connections, and particle positions:
"""
for c in self.constraints:
c.draw()
#for p in self:
# p.draw()
def lockParticle(self):
"""
If the mouse LMB is pressed for the first time on a particle, the particle
will assume the mouse motion. When it is pressed again, it will lock
the particle in space.
"""
mousePos = Vec2d(pygame.mouse.get_pos())
for p in self:
dist2mouse = sqrt(abs(pow(p.currentPos.x -
mousePos.x, 2) +
pow(p.currentPos.y -
mousePos.y, 2)))
if dist2mouse < 10:
if not p.followMouse:
p.locked = True
p.followMouse = True
p.oldPos = Vec2d(p.currentPos)
else:
p.followMouse = False
def unlockParticle(self):
"""
If the RMB is pressed on a particle, if the particle is currently
locked or being moved by the mouse, it will be 'unlocked'/stop following
the mouse.
"""
mousePos = Vec2d(pygame.mouse.get_pos())
for p in self:
dist2mouse = sqrt(abs(pow(p.currentPos.x -
mousePos.x, 2) +
pow(p.currentPos.y -
mousePos.y, 2)))
if dist2mouse < 5:
p.locked = False
#------------
# Main Program
def main():
# Screen Setup
screen = pygame.display.set_mode((WIDTH, HEIGHT))
clock = pygame.time.Clock()
# Create our grid of particles:
particleSystem = ParticleSystem(screen)
backgroundCol = Color('black')
# main loop
looping = True
while looping:
clock.tick(FRAMERATE)
pygame.display.set_caption("%s -- www.AKEric.com -- LMB: movelock - RMB: unlock - fps: %.2f"%(TITLE, clock.get_fps()) )
screen.fill(backgroundCol)
# Detect for events
for event in pygame.event.get():
if event.type == pygame.QUIT:
looping = False
elif event.type == MOUSEBUTTONDOWN:
if event.button == 1:
# See if we can make a particle follow the mouse and lock
# its position when done.
particleSystem.lockParticle()
if event.button == 3:
# Try to unlock the current particles position:
particleSystem.unlockParticle()
# Do stuff!
particleSystem.timeStep()
particleSystem.draw()
# update our display:
pygame.display.update()
#------------
# Execution from shellicon:
if __name__ == "__main__":
print "Running Python version:", sys.version
print "Running PyGame version:", pygame.ver
print "Running %s.py"%TITLE
sys.exit(main())
因为这两个程序的工作方式大致相同，但Python版本要慢得多，这让我想知道:
- 这种性能差异是Python本质的一部分吗?
- 如果我想从自己的Python程序中获得更好的性能，我应该做些什么与上面不同的事情?例如，将所有粒子的属性存储在一个数组中，而不是使用单个对象等。
编辑:回答! !
@ E先生在评论中链接的PyCon演讲，和@A。Rosa的回答和链接的资源都极大地帮助了我们更好地理解如何编写又好又快的python代码。我现在把这个页面加为书签以备将来参考:D
在Python Wiki的性能提示部分有一篇Guido van Rossum的文章链接。在它的结尾，你可以读到下面的句子:
如果你觉得需要速度，那就使用内置函数——你不能打败用c编写的循环。
文章继续列出循环优化的指导方针。我推荐这两种资源，因为它们提供了关于优化Python代码的具体和实用的建议。
在benchmarksgame.alioth.debian.org上也有一组著名的基准测试，在那里您可以找到不同机器上不同程序和语言之间的比较。可以看到，有很多变量在起作用，使得这样广泛的不可能状态Java比Python快。这通常概括为一句话"语言没有速度;实现做"。可以在代码中使用内置函数应用更python化和更快的替代方案。例如，有几个嵌套循环(其中一些不需要处理整个列表)可以使用imap或列表推导式重写。PyPy也是提高性能的另一个有趣的选项。我不是Python优化方面的专家，但是有很多非常有用的技巧(注意不要用Python写Java就是其中之一!)。
参考资料和关于SO的其他相关问题:
- Python和C的性能差异
- 是合理的集成python与c的性能?
- http://www.ibm.com/developerworks/opensource/library/os-pypy-intro/index.html?ca=drs-
- http://pyevolve.sourceforge.net/wordpress/?p=1189
如果你像编写Java那样编写Python，它当然会变慢，习惯的Java不能很好地转换为习惯的Python。
这种性能差异是Python本质的一部分吗?如果我想从我自己的Python程序中获得更好的性能，我应该做些什么与上面不同的事情?例如，将所有粒子的属性存储在一个数组中，而不是使用单个对象等。
没有看到你的代码很难说。
下面是python和java之间的差异的不完整列表，这些差异有时会影响性能:
-
Processing使用即时模式canvas，如果你想在Python中获得类似的性能，你也需要使用即时模式canvas。大多数GUI框架中的画布(包括Tkinter画布)都是保留模式，这种模式更容易使用，但固有地比即时模式慢。您需要使用像pygame, SDL或Pyglet提供的直接模式画布。
-
Python是动态语言，这意味着实例成员访问、模块成员访问和全局变量访问是在运行时解析的。在python中，实例成员访问、模块成员访问和全局变量访问实际上是字典访问。在java中，它们在编译时解决，并且从本质上讲要快得多。将频繁访问的全局变量、模块变量和属性缓存到局部变量
-
x, range()生成一个具体的列表，在python中，使用迭代器
for item in list进行迭代通常比使用迭代变量for n in range(len(list))进行迭代快。你应该几乎总是直接使用iterator进行迭代，而不是使用range(len(…))。 -
Python的数字是不可变的，这意味着任何算术计算都会分配一个新对象。这就是为什么普通python不太适合低级计算的原因之一;大多数希望能够编写低级计算而不必诉诸于编写C扩展的人通常使用cython、psyco或numpy。但是，这通常只有在进行数百万次计算时才会成为问题。
这只是部分的，非常不完整的列表，还有许多其他原因可以解释为什么将java翻译成python会产生次优代码。如果没有看到代码，就不可能知道需要做什么不同的事情。优化的python代码通常看起来与优化的java代码非常不同。
我还建议阅读其他物理引擎。有一些开源引擎使用各种方法来计算"物理"。
- Newton Game Dynamics
- 花栗鼠 子弹
- Box2D
- ODE (Open Dynamics Engine)
也有大多数引擎的端口:
- Pymunk
- PyBullet
- PyBox2D
- PyODE
如果你通读这些引擎的文档，你会经常发现声明说它们针对速度(30fps - 60fps)进行了优化。但如果你认为他们可以在计算"真正的"物理时做到这一点，那你就错了。大多数引擎计算物理到一个普通用户不能光学区分"真实"的物理行为和"模拟"的物理行为。然而，如果你调查了这个错误，如果你想编写游戏，它就可以忽略不计了。但如果你想学物理，所有这些引擎对你来说都没用。这就是为什么我会说，如果你在做一个真正的物理模拟，你会比那些引擎设计得慢，你永远不会超过另一个物理引擎。
基于粒子的物理模拟很容易转化为线性代数运算。矩阵运算。Numpy提供了这样的操作，这些操作在底层是用Fortran/C/c++实现的。编写良好的python/Numpy代码(充分利用语言&