最近在做毕业论文，涉及到光线追踪，NURBS曲面等。

我已经把NURBS曲面相关的东西弄好了。现在做光线追踪。不知道有什么好方法能求直线（光）和NURBS的交点。

我的想法是算出曲线上点和直线的距离，如果距离是0就是直线和NURBS曲线的交点了。

对于2d的NURBS曲线，只有一个参数u，我写了个程序dist（u）来算曲线上对应点到直线的距离，然后 fmin（dist，0.5）来求交点的参数u。

可是一运行就死机。。。

 

不知道坛子里有木有人做过光线追踪方面的东西。有没有比较好的算法来求直线和曲线的交点？（spline或者NURBS都行。）

把程序贴出来吧。

# -*- coding: utf-8 -*-

"""

Created on Thu Feb 16 00:50:04 2012

calculate the crosspoint a line with a curve  

 

@author: alaya

"""

 

from pylab import plot,show

from numpy import array,zeros,ones,linalg,sqrt,dot,linspace,shape,\

                  insert,append,cross,abs

 

from scipy.optimize import fmin

 

 

 

 

 

######## The NURBS ralated function ########

 

'''find span'''

def findspan(p,u,U):

       n=len(U)-p-1

       if u==U[n+1]: return (n-1);

       low = int(p) 

       high = int(n+1)

       mid = int((low + high)/2)

       while u<U[mid] or u>=U[mid+1]:

         if u<U[mid]:high = mid

         else: low = mid

         mid = int((low+high)/2)

       return (mid)

 

 

 

'''basis function'''

 

def basis(u,p,U):                 

        '''geben die nicht nur basis funktion als ein Reihe von N[0]=1 bis N[p]''' 

        i = findspan (p,u,U)

 

 

        N = array(ones(p+1))

        left = array(zeros(p+1))

        right = array(zeros(p+1))

        for  j in xrange(1,p+1):

 

           left [j] = u-U[i+1-j]

           right[j] = U[i+j]-u

           saved = 0

           for r in xrange (j):

              temp = N[r]/(right[r+1]+left[j-r])

              N[r] = saved + right[r+1]*temp

              saved = left[j-r]*temp

           N[j]= saved

 

        return (N)

 

 

def d_basis(u,p,U,k=1):  #calc devariation basisfunction first order,p72,eq 2.7

        i = findspan (p,u,U)

 

 

        N = ones([p+1,p+1])   #basisfunction

        Nd = zeros([k+1,p+1])  # derivated basisfunction, k derivation order

 

 

        left = zeros(p+1)

        right = zeros(p+1)

 

        for  j in xrange(1,p+1):

 

           left [j] = u-U[i+1-j]

           right[j] = U[i+j]-u

           saved = 0.0

           for r in xrange (j):

              N[j][r] = right[r+1]+left[j-r]

              temp = N[r][j-1]/N[j][r]

 

              N[r][j] = saved + right[r+1]*temp

              saved = left[j-r]*temp

           N[j][j]= saved

 

        Nd[0] = N[:,p]

 

        for r in xrange(p+1):

            a = ones([2,p+1]) #coefficient

            s1,s2 = 0,1

 

            for g in xrange(1,k+1):

                d=0.0

                rg = r-g

                pg = p-g

 

 

                if r>=g:

                    a[s2][0] = a[s1][0]/N[pg+1][rg]

                    d= a[s2][0]*N[rg][pg]

 

                if rg>=-1:

                    j1 = 1

                else: j1 = -rg

 

 

                if r-1<=pg:

                    j2 = g-1

                else: j2 = p-r

 

 

                for j in xrange(j1,j2+1):

                    a[s2][j] = a[s1][j]-a[s1][j]-a[s1][j-1]/N[pg+1][rg+j]

                    d += a[s2][j]*N[rg+j][pg]

                if r<=pg:

                    a[s2][g] = -a[s1][g-1]/N[pg+1][r]

                    d += a[s2][g]*N[r][pg]

                Nd[g][r] = d

                s1,s2 = s2,s1

 

        r = p

        for g in xrange(1,k+1):

            for j in xrange(p+1):     

                Nd[g][j]*=r

            r *= p-g

 

        return Nd

 

 

 

def curvepoints(p,U,P,u,W=zeros(1)):  #calculate the point on curve

        '''p: degree:::U: knot vector:::u: knot:::P:  control points:::n: the nummber of points'''

 

        span = findspan(p,u,U)

 

        N = basis(u,p,U)

 

        #print '\nN:',N,'P:',P[span-p:span+1],'\n'

 

        if len(W) != 1:

 

            Pw = P*W

            point = dot(N,Pw[span-p:span+1])

            point = point/dot(N,W[span-p:span+1])

 

        else: 

 

            point = dot(N,P[span-p:span+1])

 

        return (point)     

 

def curvderiv(p,U,P,u,g,W=zeros(3)): #calculate the dx/du on curve

        '''p: degree:::U: knot vector:::u: knot:::P:  control points

           n: the nummber of points:::W weightvector:::g derivation order'''   

        span = findspan(p,u,U)

        Nd = d_basis(u,p,U,g)

        Ndg = Nd[g]

        Nd0 = Nd[0]

        if len(W) != 3:

 

            Pw = P*W

            A =  dot (Nd0,Pw[span-p:span+1])

            derA = dot(Ndg,Pw[span-p:span+1])

            w = dot(Nd0,W[span-p:span+1])

            derw = dot(Ndg,W[span-p:span+1])

 

            der = (w*derA-derw*A)/w**2

 

        else: 

 

            der = dot(Ndg,P[span-p:span+1])

        return der

 

 

 

def add_CP_single(Pw,U,p,u):#P149

    Pw1 = zeros([shape(Pw)[0],shape(Pw)[1]+1])

    i = findspan(p,u,U) 

 

    U1 = insert(U,i+1,u)

    #print i,u,U1,'i,u,U1'   

    Pw1[:,:i-p+1] = Pw[:,:i-p+1]

    Pw1[:,i+1:] = Pw[:,i:]

    #print Pw1,'tempered Pw'

    for k in range(i-p+1,i+1):

 

        temp = (u-U[k])/(U[k+p]-U[k])

        #print k,temp

        Pw1[:,k] = temp * Pw[:,k] + (1-temp) * Pw[:,k-1]

 

    return U1,Pw1    

 

 

 

def NURBS_deg_plus(U,p): #

 

    U1 = array([])

    for i in range(len(U)-1):

 

        U1 = append(U1,U[i])

        if U[i] < U[i+1]:

            U1 = append(U1,U[i])

 

    U1 = append(U1,U[-2:])        

 

 

 

 

    return U1

 

def BEZIER_deg_plus(CPw):

    l = len(CPw[0])

    CPw1 = append(CPw,[[0],[0],[0]],axis = 1)

    for i in range(1,l):

        CPw1[:,i] = (1-float(i)/l)*CPw[:,i] + float(i)/l*CPw[:,i-1]    

    #print CPw1,'zwichen Punkt\n'

    CPw1[:,0] = CPw[:,0]

    CPw1[:,-1] = CPw[:,-1]

 

    return CPw1

 

######## END The NURBS ralated function ########

##### BEGIN The Raytrace ralated function ######

 

def dist(P0,V,P): #P0,V defined a line in 2d, P is the other point

    Vn = V/sqrt(sum(V**2))  #the distance a point to a line

    d = cross(P-P0,Vn)

    return abs(d)

 

 

'''    

def distNURBS(u):#calc the distance between a points on NURBS @ u and the line

    Px = curvepoints(p,U,CP[0],u,W)

    Py = curvepoints(p,U,CP[1],u,W)

    PN = array([Px,Py]) #points of NURBS

    d = dist(P0,V,PN)

    #d = cross(PN-P0,V)

    print d

    return abs(d)

'''        

 

def refractThrough(V,normal,n):

        Vn = V/linalg.norm(V)

        normal = normal/linalg.norm(normal)

        cosa = dot(Vn,normal)

        sina = sqrt(1-cosa**2)

        sinb = sina/n

        cosb = sqrt(1-sinb**2)

        temp = n*cosb-cosa

        Vn1 = Vn+temp*normal

        return Vn1

 

 

 

 

###### END The Raytrace ralated function #######

 

 

 

test_P0 = array([-1.0,0])

 

test_V = array([1,0.1])

test_P1 = test_P0 + 10*test_V

plot([-1,9],[0,1])

test_CP = array([[0.0,1.0,1.8,2.4,2.8,3.0],[0.0,1.0,2.0,3.0,4.0,5.0]])

plot(test_CP[0],test_CP[1])

test_U = array([0,0,0,0,0.3333,0.6667,1,1,1,1])

#Pcross,Dcross = CrossNURBS(3,test_U,test_CP,test_P0,test_V)

Px = curvepoints(3,test_U,test_CP[0],0)

Py = curvepoints(3,test_U,test_CP[1],0)

PN = array([Px,Py])

print dist(test_P0,test_V,PN),'distance u=0'

 

 

def find_u(u):

    u = u-int(u)

    tempx = curvepoints(3,test_U,test_CP[0],u)

    tempy = curvepoints(3,test_U,test_CP[1],u)

    tempP = array([tempx,tempy])

    tempd = dist(test_P0,test_V,tempP)

    return tempd**2

 

u = fmin(find_u,0.5)

print u

中间好大一段是关于NURBS曲线计算的。

我自己也在别的程序中用了fmin,都没问题。。不知道怎么这个程序一运行就死机。

稍微调试了一下，你可以试试运行find_u(-0.5)，这将使得findspan()中那个while循环无法退出。此时findspan()的参数为：

p=3,

u=-0.5

U =[ 0.      0.      0.      0.      0.3333  0.6667  1.      1.      1.    1.    ]

如果find_u的参数不可能为负数的话，你可以：

 

def find_u(u):

    if u < 0: return 100000

    ......

 谢谢楼上！还费力看代码。。

 

还想请问下，楼上是怎么调试的？

比如这中无限循环的情况系统也不会报错。。也不会告诉你哪里出问题了。。

另外有什么办法实现这种带边界限制的优化？

调试方法很简单，在find_u最前面添加一行print u，那么最后那个显示的u值就是进入无限循环的那个u值。接下来试着运行find_u(u)，程序进入死循环，按ctrl+c，可以看到强制退出程序时位置，那个就是死循环的位置。边界限制就通过在find_u中添加判断边界，如果出界则返回一个很大的数即可。

哦～～！

另外怎么才能调用fmin的时候不输出优化过程的信息比如

Optimization terminated successfully.

         Current function value: 0.000000

         Iterations: 12

         Function evaluations: 24

 

 

问题已经解决了。。再次感谢。