program blur
  use pgmio

  implicit none

  double precision, parameter :: t=10.0d0, deltaT=0.2d0, k=5.0d0
  character*(*), parameter :: input='dif_tomography.pgm', output='dif_tomography_perona.pgm'

  double precision, allocatable :: u(:,:), nu(:,:)
  double precision :: north, south, east, west, level
  integer :: w, h, offset, n, i, j

  call pgmsize(input, w, h, offset)

  allocate(u(0:w+1,0:h+1))
  u=0

  allocate(nu(w,h))

  call pgmread(input, offset, w, h, u, 0, 0)

  level = 0
  do while (level<t)
    do j=1,h
        do i=1,w
            north = u(i-1,j)-u(i,j)
            south = u(i+1,j)-u(i,j)
            east = u(i,j+1)-u(i,j)
            west = u(i,j-1)-u(i,j)
            nu(i,j) = u(i,j) + deltaT*(g(north)*north+g(south)*south+g(east)*east+g(west)*west)
        end do
    end do
    u(1:w,1:h) = nu(1:w,1:h)
    level = level + deltaT
  end do

  deallocate(u)

  call pgmwriteheader(output, w, h)
  call pgmappendbytes(output, nu, 1, 1)

  deallocate(nu)

  contains
      function g(x) result(v)
        implicit none
        double precision, intent(in) :: x
        double precision :: v
        v = 1/(1+(x/k)**2)
      end function g

end program blur