computeCcorr :: ImageGS -> ImageGS -> HIP.Image HIP.RPU HIP.Y Double
computeCcorr img tmpl = HIP.makeImage resultDims computePx
  where
    computePx :: (Int, Int) -> HIP.Pixel HIP.Y Double
    computePx coord = HIP.PixelY $ computePxNumer coord / sqrt (sumTmpl * sumImg)
      where
        sumTmpl = computeSqSum tmpl (0,0)
        sumImg = computeSqSum img coord
    computePxNumer :: (Int, Int) -> Double
    computePxNumer (dRow, dCol) =
      sum $ do
          r <- [0 .. tmplRows - 1]
          c <- [0 .. tmplCols - 1]
          let imgPx = HIP.getPxC (HIP.index img (dRow + r, dCol + c)) HIP.LumaY
          let tmplPx = HIP.getPxC (HIP.index tmpl (r, c)) HIP.LumaY
          pure $ imgPx * tmplPx
    resultDims = (imgRows - tmplRows + 1, imgCols - tmplCols + 1)
    (imgRows, imgCols) = HIP.dims img
    (tmplRows, tmplCols) = HIP.dims tmpl
    computeSqSum imgX (offRow, offCol) = sum $ do
      r <- [0 .. tmplRows - 1]
      c <- [0 .. tmplCols - 1]
      let v = HIP.getPxC (HIP.index imgX (offRow + r, offCol + c)) HIP.LumaY
      pure $ v * v
computeCcorr :: ImageGS -> ImageGS -> HIP.Image HIP.RPU HIP.Y Double
computeCcorr img tmpl = HIP.makeImage resultDims computePx
  where
    computePx :: (Int, Int) -> HIP.Pixel HIP.Y Double
    computePx coord = HIP.PixelY $ computePxNumer coord / sqrt (sumTmpl * sumImg)
      where
        sumTmpl = computeSqSum tmpl (0,0)
        sumImg = computeSqSum img coord
    computePxNumer :: (Int, Int) -> Double
    computePxNumer (dRow, dCol) =
      sum $ do
          r <- [0 .. tmplRows - 1]
          c <- [0 .. tmplCols - 1]
          let imgPx = HIP.getPxC (HIP.index img (dRow + r, dCol + c)) HIP.LumaY
          let tmplPx = HIP.getPxC (HIP.index tmpl (r, c)) HIP.LumaY
          pure $ imgPx * tmplPx
    resultDims = (imgRows - tmplRows + 1, imgCols - tmplCols + 1)
    (imgRows, imgCols) = HIP.dims img
    (tmplRows, tmplCols) = HIP.dims tmpl
    computeSqSum imgX (offRow, offCol) = sum $ do
      r <- [0 .. tmplRows - 1]
      c <- [0 .. tmplCols - 1]
      let v = HIP.getPxC (HIP.index imgX (offRow + r, offCol + c)) HIP.LumaY
      pure $ v * v
computeCcorr img tmpl = HIP.makeImage resultDims computePx
  where
    computePx :: (Int, Int) -> HIP.Pixel HIP.Y Double
    computePx coord = HIP.PixelY $ computePxNumer coord / sqrt (sumTmpl * sumImg)
      where
        sumTmpl = computeSqSum tmpl (0,0)
        sumImg = computeSqSum img coord
    computePxNumer :: (Int, Int) -> Double
    computePxNumer (dRow, dCol) =
      sum $ do
          r <- [0 .. tmplRows - 1]
          c <- [0 .. tmplCols - 1]
          let imgPx = HIP.getPxC (HIP.index img (dRow + r, dCol + c)) HIP.LumaY
          let tmplPx = HIP.getPxC (HIP.index tmpl (r, c)) HIP.LumaY
          pure $ imgPx * tmplPx
    resultDims = (imgRows - tmplRows + 1, imgCols - tmplCols + 1)
    (imgRows, imgCols) = HIP.dims img
    (tmplRows, tmplCols) = HIP.dims tmpl
    computeSqSum imgX (offRow, offCol) = sum $ do
      r <- [0 .. tmplRows - 1]
      c <- [0 .. tmplCols - 1]
      let v = HIP.getPxC (HIP.index imgX (offRow + r, offCol + c)) HIP.LumaY
      pure $ v * v