2013-05-09 96 views
0

我寫了一些代碼,使我能夠從DICOM文件中獲取數據並將數據分離爲單獨的FID信號。「參數必須是可調用函數」當使用scipy.integrate.quad

from pylab import * 
import dicom 
import numpy as np 
import scipy as sp 

plan=dicom.read_file("1.3.46.670589.11.38085.5.22.3.1.4792.2013050818105496124") 
all_points = array(plan.SpectroscopyData) 
cmplx_data = all_points[0::2] -1j*all_points[1::2] 
frames = int(plan.NumberOfFrames) 
fid_pts = len(cmplx_data)/frames 
del_t = plan.AcquisitionDuration/(frames * fid_pts) 

fid_list = [] 
for fidN in arange(frames): 
    offset = fidN * fid_pts 
    current_fid = cmplx_data[offset:offset+fid_pts] 
    fid_list.append(current_fid) 

我現在想確定一些這個數據的話,應用傅里葉變換和移位我試着使用SciPy的的四功能,遇到了以下錯誤後:

spec = fftshift(fft(fid_list[0]))) 
sp.integrate.quad(spec, 660.0, 700.0) 



error          Traceback (most recent call last) 
/home/dominicc/Experiments/In Vitro/Glu, Cr/Phantom 1: 10mM Cr, 5mM Glu/WIP_SV_PRESS_ME_128TEs_5ms_spacing_1828/<ipython-input-107-17cb50e45927> in <module>() 
----> 1 sp.integrate.quad(fid, 660.0, 700.0) 

/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst) 
243  if type(args) != type(()): args = (args,) 
244  if (weight is None): 
--> 245   retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points) 
246  else: 
247   retval = _quad_weight(func,a,b,args,full_output,epsabs,epsrel,limlst,limit,maxp1,weight,wvar,wopts) 

/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.pyc in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points) 
307  if points is None: 
308   if infbounds == 0: 
--> 309    return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit) 
310   else: 
311    return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit) 

error: First argument must be a callable function. 
莫非

有人請建議一種方法來使這個對象可調用?閱讀What is a "callable" in Python?後,我還是不太清楚。讀完後,我試了

def fid(): 
    spec = fftshift(fft(fid_list[0])) 
    return spec 

哪個返回了相同的錯誤。

任何幫助將不勝感激,謝謝。

+0

它不應該是一個需要參數的函數嗎?還有什麼可以整合? – 2013-05-09 13:42:45

+0

這個問題沒有意義。你不能使用'scipy.integrate.quad'集成一個數組。它集成了一個*函數*,它接受一個參數(例如'sin(x)')。使用類似'scipy.integrate.simps'的東西來集成一組函數值。 – talonmies 2013-05-09 13:44:11

回答

2

由於您正在整合僅在網格上定義的函數(即數字數組),因此您需要使用「Integration, given fixed samples」中的一個例程。

事實上,第一個參數quad()must be a function,你可以調用(「可調用」)的東西。舉例來說,你可以這樣做:

>>> from scipy.integrate import quad 
>>> def f(x): 
...  return x**2 
... 
>>> quad(f, 0, 1) 
(0.33333333333333337, 3.700743415417189e-15) 

這是可以做到的,而不是抽樣f()

>>> from scipy.integrate import simps 
>>> x_grid = numpy.linspace(0, 1, 101) # 101 numbers between 0 and 1 
>>> simps(x_grid**2, dx=x_grid[1]-x_grid[0]) # x_grid**2 is an *array*. dx=0.01 between x_grid values 
0.33333333333333337 

你的情況就像是在第二個例子:你整合採樣功能,因此很自然地使用一個或cumtrapz(),simps()romb()

相關問題