ImageConvolve.scala

// import de.sciss.file._

import edu.emory.mathcs.jtransforms.fft.DoubleFFT_2D
import javax.imageio.ImageIO
import java.awt.image.BufferedImage

val dir = userHome/"Documents"/"projects"/"Anemone"/"ConvTest"
val p1  = dir/"Rattle_shot150131sz.png"
val p2  = dir/"Wolkenpumpe141214_135413sz.png"

def extractChannel(in: BufferedImage, chan: Int): Array[Array[Double]] = {
  val shift = chan * 8
  Array.tabulate(in.getHeight) { y => 
    Array.tabulate(in.getWidth) { x =>
      val i = (in.getRGB(x, y) >>> shift) & 0xFF
      i.toDouble / 0xFF
    }
  }
}

def fillChannel(in: Array[Array[Double]], out: BufferedImage, chan: Int): Unit = {
  val shift = chan * 8
  val mask  = ~(0xFF << shift)
  for (y <- 0 until in.length) {
    val v = in(y)
    for (x <- 0 until v.length) {
      val d = v(x)
      val i = (d.clip(0, 1) * 0xFF + 0.5).toInt << shift
      val j = out.getRGB(x, y)
      val k = j & mask | i
      out.setRGB(x, y, k)
    }
  }
}

def mulC(a: Array[Array[Double]], b: Array[Array[Double]], scale: Double = 1.0): Unit = {
  for (y <- 0 until a.length) {
    val va = a(y)
    val vb = b(y)
    for (x <- 0 until va.length by 2) {
      val re = va(x) * vb(x)   - va(x+1) * vb(x+1)
      val im = va(x) * vb(x+1) + va(x+1) * vb(x)
      va(x)   = re * scale
      va(x+1) = im * scale
    }
  }
}

val pOut = dir / "out.png"

def perform(div: Double = 256): Unit = {
  val i1  = ImageIO.read(p1)
  val i2  = ImageIO.read(p2)

  val w = i1.getWidth
  val h = i1.getHeight
  require(w == i2.getWidth && h == i2.getHeight)
  require(w.isPowerOfTwo && h.isPowerOfTwo)
  
  val fft = new DoubleFFT_2D(h, w)

  for (chan <- 0 until 3) {
    val b1 = extractChannel(i1, chan)
    val b2 = extractChannel(i2, chan)
    fft.realForward(b1)
    fft.realForward(b2)
    mulC(b1, b2, scale = 1.0/div)
    fft.realInverse(b1, true)
    fillChannel(b1, i1, chan = chan)
  }
  ImageIO.write(i1, "png", pOut)
}

perform(2048*3/2)