周波数フィルタリング

フーリエ変換

ローパスフィルタ・ハイパスフィルタ・バンドパスフィルタ

library( imager )

clamp = function( array, min, max ){
  array[ array < min ] = min
  array[ array > max ] = max
  return( array )
}

swapPixels = function( img, reverse = FALSE ){
  a = as.array( img )
  w = width( img )
  h = height( img )
  midx.minus1 = floor( w / 2 )
  midy.minus1 = floor( h / 2 )
  if( mod( w, 2 ) == 0 ){
    left = a[ 1:midx.minus1+0, , , ]
    a[ 1:midx.minus1+0, , , ] = a[ (midx.minus1+1):w, , ,  ]
    a[ (midx.minus1 + 1):w, , , ] = left
  } else {
    if( reverse ){
      left = a[ 1:midx.minus1+0, , , ]
      a[ 1:(midx.minus1+1), , , ] = a[ (midx.minus1+1):w, , ,  ]
      a[ (midx.minus1+2):w, , , ] = left
    } else {
      left = a[ 1:(midx.minus1+1), , , ]
      a[ 1:midx.minus1+0, , , ] = a[ (midx.minus1+2):w, , ,  ]
      a[ (midx.minus1+1):w, , , ] = left
    }
  }
  if( mod( h, 2 ) == 0 ){
    top = a[ , 1:midy.minus1, , ]
    a[ , 1:midy.minus1, , ] = a[ , (midy.minus1 + 1):h, , ]
    a[ , (midy.minus1 + 1):h, , ] = top
  } else {
    if( reverse ){
      top = a[ , 1:midy.minus1, , ]
      a[ , 1:(midy.minus1 + 1), , ] = a[ , (midy.minus1 + 1):h, , ]
      a[ , (midy.minus1 + 2):h, , ] = top
    } else {
      top = a[ , 1:(midy.minus1 + 1), , ]
      a[ , 1:(midy.minus1 + 0), , ] = a[ , (midy.minus1 + 2):h, , ]
      a[ , (midy.minus1 + 1):h, , ] = top
    }
  }
  img = as.cimg( a, dim = dim( img ) )
  return( img )
}

highpass = function( img, cyclePerImage, smoothing = F ){
  img = swapPixels( img, reverse = F )
  a = as.vector( img )
  w = width( img )
  h = height( img )
  origin.x = floor( w / 2 ) + 1
  origin.y = floor( h / 2 ) + 1
  radius2 = cyclePerImage^2
  for( x in 1:w ){
    for( y in 1:h ){
      x2 = x - w / 2 - 1
      y2 = y - h / 2 - 1
      if( smoothing ){
        if( !(x == origin.x & y == origin.y) ){
          sigma = 1/(2*pi*cyclePerImage)
          weight = 1 - exp(-2*pi*sigma^2*(x2^2+y2^2))
          a[ x + (y-1) * w ] = weight * a[ x + (y-1) * w ]
        }
      } else {
        if( x2^2 / radius2 + y2^2 / radius2 < 1 ){
          if( !(x == origin.x & y == origin.y) ){
            a[ x + (y-1) * w ] = 0
          }
        }
      }
    }
  }
  img = as.cimg( a, x = w, y = h, z = 1, cc = 1 )
  img = swapPixels( img, reverse = T )
  return( img )
}

lowpass = function( img, cyclePerImage, smoothing = F ){
  img = swapPixels( img, reverse = F )
  a = as.vector( img )
  w = width( img )
  h = height( img )
  origin.x = floor( w / 2 ) + 1
  origin.y = floor( h / 2 ) + 1
  radius2 = cyclePerImage^2
  for( x in 1:w ){
    for( y in 1:h ){
      x2 = x - w / 2 - 1
      y2 = y - h / 2 - 1
      if( smoothing ){
        if( !(x == origin.x & y == origin.y) ){
          sigma = 1/(2*pi*cyclePerImage)
          weight = exp(-2*pi*sigma^2*(x2^2+y2^2))
          a[ x + (y-1) * w ] = weight * a[ x + (y-1) * w ]
        }
      } else {
        if( x2^2 / radius2 + y2^2 / radius2 > 1 ){
          a[ x + (y-1) * w ] = 0
        }
      }
    }
  }
  img = as.cimg( a, x = w, y = h, z = 1, cc = 1 )
  img = swapPixels( img, reverse = T )
  return( img )
}

bandpass = function( ff, cyclePerImage, bandwidth = 2, smoothing = F ){
  low = cyclePerImage * 2^(-bandwidth/2)
  high = cyclePerImage * 2^(+bandwidth/2)
  ff = lapply( ff, highpass, cyclePerImage = low, smoothing = smoothing )
  ff = lapply( ff, lowpass, cyclePerImage = high, smoothing = smoothing )
  return( ff )
}
as.fftMagnitude = function( ff, clamp = F ){
  p = sqrt( ff$real^2 + ff$imag^2 )
  if( clamp ){
    max_ = min( quantile( p, (10^(1:9)-1)/10^(1:9) )[ quantile( p, (10^(1:9)-1)/10^(1:9) ) > 0 ] )
    p = clamp( p, min = 0, max = max_ )
    return( p )
  } else {
    return( p )
  }
}

as.fftPhase = function( ff, clamp = F ){
  p = atan2( ff$imag, ff$real )
  if( clamp ){
    max_ = min( quantile( p, (10^(1:9)-1)/10^(1:9) )[ quantile( p, (10^(1:9)-1)/10^(1:9) ) > 0 ] )
    p = clamp( p, min = 0, max = max_ )
    return( p )
  } else {
    return( p )
  }
}

as.phaseScrampled = function( img ){
  ff = FFT( img )
  ff.magnitude = as.fftMagnitude( ff )
  # ff.phase = as.fftPhase( ff )
  ff.phase2 = as.cimg( (runif(width(img)*height(img))-.5)/.5*pi, x = width(img), y = height(img), z = 1, cc = 1 )
  ff.real = ff.magnitude * cos( ff.phase2 )
  ff.imag = ff.magnitude * sin( ff.phase2 )
  check = FFT( ff.real, ff.imag, inverse = TRUE )$real
  return( check )
  # m0 = as.fftMagnitude( ff )
  # m1 = as.fftMagnitude( FFT( check ) )
  # m = sqrt( ff.real^2 + ff.imag^2 )
  # print( mean( ( m0 - m1 )^2 ) )
  # print( mean( ( m0 - m )^2 ) )
  # mean( ( check - img )^2 )
  # layout(t(1:4))
  # plot( img, rescale = F )
  # plot( check, rescale = T )
  # plot( swapPixels( as.fftMagnitude( ff, clamp = T ) ), rescale = T, main="Power spectrum", cex.main = 2 )
  # plot( swapPixels( as.fftMagnitude( FFT( check ), clamp = T ) ), rescale = T, main="Power spectrum", cex.main = 2 )
  # plot( img, rescale = F )
}

freqFilter = function( img, scheme, cyclePerImage, smoothing = F, bandwidth = 2 ){
  if( scheme == "highpass" ){
    ff = lapply( FFT( img ), highpass, cyclePerImage = cyclePerImage, smoothing = smoothing )
  } else if( scheme == "lowpass" ){
    ff = lapply( FFT( img ), lowpass, cyclePerImage = cyclePerImage, smoothing = smoothing )
  } else if( scheme == "bandpass" ){
    ff = bandpass( FFT( img ), cyclePerImage = cyclePerImage, smoothing = smoothing, bandwidth = bandwidth )
  } else {
    return( img )
  }
  out = FFT( ff$real, ff$imag, inverse = TRUE )$real
  return( out )
}

#
inputImg = boats
img = grayscale( inputImg )
low = freqFilter( img, scheme = "lowpass", cyclePerImage = 6, smoothing = T )
high = freqFilter( img, scheme = "highpass", cyclePerImage = 28, smoothing = T )
band = freqFilter( img, scheme = "bandpass", cyclePerImage = 15, smoothing = T, bandwidth = 2 )
plot( imlist( input = inputImg, lowpass = low, highpass = high, bandpass = band ) )