Fortran to Matlab, B-Spline,

Question

STamer 2013 年 7 月 11 日

0
リンク

この質問への直接リンク

https://jp.mathworks.com/matlabcentral/answers/81741-fortran-to-matlab-b-spline

program bspline_fit

    use toolbox
    implicit none
    integer :: ord,typ,ni,no
    real :: t,dt
    real, allocatable :: nvo(:),nvi(:),crd0(:,:),crd1(:,:)
    integer :: i,j
    ! input data (the number of control points will be equal to the number of data points)
    write (*,*) "Enter the number of data points"
    read (*,*) ni
    write (*,*) "Enter the order of the curve"
    read (*,*) ord
    no = ni+ord
    write (*,*) "Enter the type of the knot vector"
    write (*,*) "(1) periodic uniform (2) open uniform (3) nonuniform"
    read (*,*) typ
    allocate (crd0(ni,2),crd1(ni,2),nvi(ni),nvo(no))
    call data_read(ni,crd0) ! input data
    call data_vec(ni,crd0,nvi) ! parametric coordinates of the data points
    call knot_vec(typ,ord,ni,no,nvo) ! knot vector
    call control_vec(ord,ni,no,nvi,nvo,crd0,crd1) ! parametric coordinates of the control points
    call curve(ord,ni,no,nvi,nvo,crd1) ! resulting curve

end program bspline_fit

--------------------------------------

module toolbox

    contains
    subroutine data_read(n,crd)
    implicit none
    integer,intent(in) :: n
    real,intent(out) :: crd(n,2)
    integer :: i,j,k
    open (1,file = 'in.dat')
    do i=1,n ; read (1,*) k,(crd(i,j),j=1,2) ; end do
    end subroutine data_read
    subroutine data_vec(n,crd,nvi)
    implicit none
    integer,intent(in) :: n
    real,intent(in) :: crd(n,2)
    real,intent(out) :: nvi(n)
    real :: dif(2),clen(n-1)
    integer :: i
    do i=1,n-1
        dif = crd(i+1,:)-crd(i,:)
        clen(i) = sqrt(dif(1)**2+dif(2)**2)
    end do
    nvi(1) = 0
    do i=2,n
        nvi(i) = nvi(i-1)+clen(i-1)
    end do
    nvi = nvi/sum(clen)
    write (*,'(20(f5.2,2x))') (nvi(i),i=1,n) ; write (*,*)
    end subroutine data_vec
    subroutine knot_vec(typ,ord,ni,no,nvo)
    integer,intent(in) :: typ,ord,ni,no
    real,intent(out) :: nvo(no)
    integer :: i
    select case(typ)
    case(1)
    case(2)
        do i=1,no
            if (i <= ord) then
                nvo(i) = 0
            else if (i >= ni+1) then
                nvo(i) = ni+1-ord
            else
                nvo(i) = i-ord
            end if
        end do
    case(3)
    end select
    nvo = nvo/nvo(no)
    write (*,'(20(f5.2,2x))') (nvo(i),i=1,no)
    end subroutine knot_vec
    subroutine control_vec(ord,ni,no,nvi,nvo,crd0,crd1)
    integer,intent(in) :: ord,ni,no
    real,intent(in) :: nvi(ni),nvo(no),crd0(ni,2)
    real,intent(out) :: crd1(ni,2)
    real :: cmat(ni,ni)
    integer :: i
    open (3,file='tmp.dat')
    do i=1,ni
        call basis_fnc(ord,ni,no,nvo,nvi(i),cmat(i,:))
    end do
    do i=1,ni ; write (*,'(20(f6.3,2x))') (cmat(i,j),j=1,ni) ; end do ; write (3,*)
    call mat_inv(ni,cmat) ; crd1 = matmul(cmat,crd0)
    do i=1,ni ; write (3,'(i3,1x,3(f6.3,2x))') i,crd1(i,1),crd1(i,2) ; end do ; write (3,*)
    end subroutine control_vec
    subroutine basis_fnc(ord,ni,no,x,t,cv)
    implicit none
    integer,intent(in) :: ord,ni,no
    real,intent(in) :: x(no),t
    real,intent(out) :: cv(ni)
    real :: bf(ord,no-1),r
    integer :: i,j,k
    bf = 0
    k = no-1
    do i=1,ord
        if (i == 1) then
            do j=1,k
                if (t >= x(j) .and. t <= x(j+1)) then
                    bf(i,j) = 1
                else
                    bf(i,j) = 0
                end if
            end do
        else
            k = k-1
            do j=1,k
                r = x(i+j-1)-x(j) ; if (r /= 0) bf(i,j) = (t-x(j))*bf(i-1,j)/r
                r = x(i+j)-x(j+1) ; if (r /= 0) bf(i,j) = bf(i,j)+(x(i+j)-t)*bf(i-1,j+1)/r
            end do
        end if
        write (3,'(20(f6.3,2x))') (bf(i,j),j=1,k)
    end do
    cv = bf(ord,1:ni)
    end subroutine basis_fnc
    subroutine curve(ord,ni,no,nvi,nvo,crd1)
    implicit none
    integer,intent(in) :: ord,ni,no
    real,intent(in) :: nvi(ni),nvo(no),crd1(ni,2)
    real :: cv(ni),crd(2),t
    integer :: i,j,k
    real,parameter :: dt = 1./20
    open (2,file = 'out.dat')
    t = 0 ; i = 1
    do
        call basis_fnc(ord,ni,no,nvo,t,cv)
        do j=1,2 ; crd(j) = dot_product(cv,crd1(:,j)) ; end do
        write (2,'(i3,2x,2(f5.2,2x))') i,(crd(j),j=1,2)
        i = i+1 ; t = t+dt ; if (t > nvi(ni)) exit
    end do
    end subroutine curve
    ! matrix inverse
    subroutine mat_inv(n,a)
    implicit none
    integer,intent(in) :: n
    real,intent(inout) :: a(n,n)
    integer :: i
    real :: v1(n),v2(n)
    call sgetrf(n,n,a,n,v1,i) ! lapack lu decomposition
    call sgetri(n,a,n,v1,v2,n,i) ! lapack matrix inverse
    end subroutine mat_inv