nldifffilter.f90

这是一个SUSAN算法实现, 对图像处理中的边,角等特征保留和图像降噪方面有着重要的应用

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! mex file for nlDiffFilter
!
! B = nlDiffFilter(A,table,mask)
!
! INPUT(all inputs should be of type double):
!       A (m,n) - the input image
!       table(2*L+1) - a lookup table for a nonlinear function f(x)
!         width -L < x < L where L is the largest possible 
!         difference in intensities for a given an image type
!       mask(w1,w2) - a mask
!       interp (optional) - pass a fourth argument 'linear' for 
!         linear intperolation between the values in the lookup table 
!         default is 'nearest' interpolation
!
! OUTPUT:
!       B - the filtered image
!
! A general purpose non-linear filter based on the 
! differences in intensities in a neighborhood.
!
! The form of the filter is (in psuedo code):
!
! B_ij = mask_kl*f(A_ij-windowA_kl)
!
! where for every i,j one sums over all window indexes k,l
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine mexFunction(nlhs, plhs, nrhs, prhs)
 use mexf90                    ! API function definitions

 implicit none

 integer, intent(in) :: nlhs, nrhs
 integer, intent(in), dimension(*) :: prhs
 integer, intent(out), dimension(*) :: plhs

 integer :: m,n,w1,w2,L
 integer, pointer :: A, B, mask, table
 logical :: interp = .true.

 ! Check input arguments
 if(nrhs /= 4) then
    interp = .false.
    if(nrhs /= 3) then
        call mexErrMsgTxt('Function requires four input arguments.');
    end if
 end if

 ! Get data and size of the input matrix
 A => mxGetPr(prhs(1))
 m = mxGetM(prhs(1))
 n = mxGetN(prhs(1))

 ! Get the data in the lookup table
 table => mxGetPr(prhs(2))
 L = mxGetN(prhs(2))

 !get the mask data
 mask => mxGetPr(prhs(3))
 w1 = mxGetM(prhs(3))
 w2 = mxGetN(prhs(3))

 ! Create output matrix
 plhs(1) = mxCreateDoubleMatrix(m,n,0)
 B => mxGetPr(plhs(1))

 ! Call subroutine for multiplication
 call nlDiffFilt(A,B,table,mask,m,n,w1,w2,L,interp)

end subroutine mexFunction

subroutine nlDiffFilt(A,B,table,mask,m,n,w1,w2,L,interp)

  implicit none

  integer :: m,n,w1,w2,L,i,j,w1w2,r1,r2,lb,ub
  ! Now define the matrices with the actual data type and dimension
  double precision, dimension(m+(w1-1),n+(w2-1)) :: tempA
  double precision, dimension(m,n) :: A,B
  double precision, dimension(w1,w2) :: mask
  double precision, dimension(-(L-1)/2:(L-1)/2) :: table
  integer, dimension(w1*w2) :: temp1
  double precision, dimension(w1*w2) :: temp2,f,vMask
  logical :: interp
  
  w1w2 = w1*w2
  r1 = (w1-1)/2
  r2 = (w2-1)/2

  vMask = reshape(mask,(/w1*w2/))

  tempA = reshape(spread(0.,1,m*n+2*m*r2+2*n*r1+4*r1*r2),(/m+2*r2,n+2*r2/))
  tempA(1+r1:m+1+r1,1+r2:1+n+r2) = A

  if(.not.interp) then

  do i = 1+r1,r1+m
       do j = 1+r2,r2+n
           !this horrendous operation deserves an explantion
           !  1) access a subset of tempB corresponding to the window
           !  2) subtract the entire subsection from tempB(i,j)
           !  3) make result into integral indexes and reshape 
           !  4) lookup function value through indexes and assign
           !  5) mutiple this vector by the mask and take the sum and assign to B

           temp1 = reshape(int(tempA(i-r1:i+r1,j-r2:j+r2) - tempA(i,j)),(/w1w2/))
           B(i-r1,j-r2) = sum(vMask*table(temp1))
       end do
  end do
  else

  do i = 1+r1,r1+m
       do j = 1+r2,r2+n
       
           !this horrendous operation deserves an explantion
           !  1) access a subset of tempB corresponding to the window
           !  2) subtract the entire subsection from tempB(i,j)
           !  3) make result into integral indexes and reshape 
           !  4) lookup function value through indexes and assign
           !  5) mutiple this vector by the mask and take the sum and assign to B

           temp2 = reshape(tempA(i-r1:i+r1,j-r2:j+r2) - tempA(i,j),(/w1w2/))
           temp1 = ceiling(temp2+tiny(1.))
           f = temp1-temp2
 
           B(i-r1,j-r2) = sum(vMask*((1-f)*table(temp1)+f*table(temp1-1)))
       end do
  end do
  
  end if

end subroutine nlDiffFilt