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 xTechNotes.blogspot.com
! *
! *   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
! *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! *    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)
                    endif                      
                    enddo
                enddo
            enddo
        enddo
        return
    end subroutine conv2