Friday, November 12, 2010

Matlab Function in Fortran - conv2 - Convolution in 2D

Below is the Fortran implementation of the Convolution 2D algorithm, or conv2 as found in Matlab. There are many useful utility functions in Matlab and occassionally an important but rare (difficult to find source in the internet) algorithm such as this convolution algorithm. Hopefully there will be more Fortran version of Matlab utility functions appearing here in the future.

The one listed below is unpolished and unoptimized and not suitable for large inputs. It is a very basic implementation taken directly from the mathematical definition of discrete convolution in 2D. It is also tested only for a small number of cases. Corrections, suggestions and comments are welcomed.

! * CONV2 - Convolution in 2D
! *
! *     (C)  Copyright 2010
! *
! *   This program is free software: you can redistribute it and/or modify
! *   it under the terms of the GNU General Public License as published by
! *   the Free Software Foundation, either version 3 of the License, or
! *   (at your option) any later version.
! *
! *    This program is distributed in the hope that it will be useful,
! *    but WITHOUT ANY WARRANTY; without even the implied warranty of
! *    GNU General Public License for more details.
! *
! *    You should have received a copy of the GNU General Public License
! *    along with this program.  If not, see .
! */  
    subroutine conv2(nXi, nXj, nHi, nHj, nYi, nYj, X, H, Y)
        INTEGER, INTENT(IN) :: nXi, nXj, nHi, nHj           ! Dimensions of Input, Kernel
        INTEGER, INTENT(IN) :: nYi, nYj                     ! Dimensions of Output - MUST BE nXi+nHi-1, nXj+nHj-1
        REAL, INTENT(IN) :: X(0:nXi-1, 0:nXj-1)   ! Input
        REAL, INTENT(IN) :: H(0:nHi-1, 0:nHj-1)   ! Kernel
        REAL, INTENT(OUT) :: Y(0:nYi-1, 0:nYj-1)  ! Output
        integer :: im, in, ii, ij, ip, iq
        Y = 0.0d0
        if( nYi .ne. nXi+nHi-1  .or. nYj .ne. nXj+nHj-1 ) RETURN
        do in = 0, nYj-1
            do im = 0, nYi-1
                do ij = 0, nXj-1
                    do ii = 0, nXi-1
                    ip = im - ii
                    iq = in - ij
                    if (ip .ge. 0 .and. ip .le. nHi-1 .and. iq .ge. 0 .and. iq .le. nHj-1) then
                        Y(im, in) = Y(im, in) + X(ii, ij) * H(ip, iq)
    end subroutine conv2

No comments: