!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 mexFunctionsubroutine 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 ifend subroutine nlDiffFilt