Improve perofrmance of interpolation by inlining

hip/bench/Resize.hs

@@ -7,7 +7,6 @@ import Codec.Picture.Extra
 import Criterion.Main
 import qualified Data.Massiv.Array.IO as M
 import Graphics.Image as I
-import Graphics.Image.Processing.Filter
 import Prelude as P
 
 main :: IO ()
@@ -16,16 +15,31 @@ main = do
       up = Sz (125 :. 400)
   defaultMain
     [ env (readImageRGB8 "images/frog.jpg") $ \img ->
-        bgroup "Resize" [benchScale down img, benchScale up img]
+        bgroup
+          "Resize"
+          [ benchBilinearResize down img
+          , benchBilinearResize up img
+          , benchBicubicResize down img
+          , benchBicubicResize up img
+          ]
     ]
 
 
-benchScale :: Sz2 -> Image (SRGB 'Linear) Word8 -> Benchmark
-benchScale sz@(Sz2 m n) ~img@(Image a) =
+benchBilinearResize :: Sz2 -> Image (SRGB 'Linear) Word8 -> Benchmark
+benchBilinearResize sz@(Sz2 m n) ~img@(Image a) =
   bgroup
     ("Bilinear " ++ show sz)
-    [ bench "HIP" $ nf (resize Bilinear (Fill 0) sz) img
+    [ bench "HIP" $ nf (resize Bilinear (Fill 0) sz) i
     , bench "JuicyPixels-extra" $ nf (scaleBilinear n m) jp
     ]
   where
+    i = I.map (fmap toFloat) img
     jp = M.toJPImageRGB8 (M.toImageBaseModel a)
+
+benchBicubicResize :: Sz2 -> Image (SRGB 'Linear) Word8 -> Benchmark
+benchBicubicResize sz img =
+  bgroup
+    ("Bcubic " ++ show sz)
+    [bench "HIP" $ nf (resize (Bicubic 0.5) (Fill 0) sz) i]
+  where
+    i = I.map (fmap toFloat) img

hip/src/Graphics/Image/Processing/Geometric.hs

@@ -294,10 +294,10 @@ rotate270 = transpose . flipH
 --
 -- <<images/frog.jpg>> <<images/frog_rotate330.png>>
 --
-rotate :: (ColorModel cs e, Interpolation method) =>
+rotate :: (RealFloat e, ColorModel cs e, Interpolation method) =>
           method -- ^ Interpolation method to be used
        -> Border (Pixel cs e) -- ^ Border handling strategy
-       -> Double -- ^ Angle Θ in radians
+       -> e -- ^ Angle Θ in radians
        -> Image cs e -- ^ Source image
        -> Image cs e -- ^ Rotated image
 rotate method border theta' (Image arr) =
@@ -307,7 +307,7 @@ rotate method border theta' (Image arr) =
         !j' = jD * cosTheta - iD * sinTheta - 0.5
     in interpolate method (A.handleBorderIndex border sz (A.index' arr)) (i', j')
   where
-    !theta = angle0to2pi (-theta') -- invert angle direction and put it into [0, 2*pi) range
+    !theta = angle0to2pi (- theta') -- invert angle direction and put it into [0, 2*pi) range
     !sz@(Sz2 m n) = A.size arr
     !mD = fromIntegral m
     !nD = fromIntegral n
@@ -334,59 +334,58 @@ rotate method border theta' (Image arr) =
 --
 -- <<images/frog.jpg>> <<images/frog_resize.jpg>>
 --
-resize :: (ColorModel cs e, Interpolation method) =>
+resize :: (RealFloat e, ColorModel cs e, Interpolation method) =>
           method -- ^ Interpolation method to be used during scaling.
        -> Border (Pixel cs e) -- ^ Border handling strategy
        -> Sz2   -- ^ Dimensions of a result image.
        -> Image cs e -- ^ Source image.
        -> Image cs e -- ^ Result image.
-resize method border sz'@(Sz2 m' n') (Image arr) = --Image $ A.compute warr
-  Image (A.makeArray (A.getComp arr) sz' getNewPx)
-  where
-    sz@(Sz2 m n) = A.size arr
-    !fM = fromIntegral m' / fromIntegral m
-    !fN = fromIntegral n' / fromIntegral n
-    getNewPx (i :. j) =
-      interpolate
-        method
-        (A.handleBorderIndex border sz (A.index' arr))
-        ( (fromIntegral i + 0.5) / fM - 0.5
-        , (fromIntegral j + 0.5) / fN - (0.5 :: Double))
-    {-# INLINE getNewPx #-}
+resize method border sz'@(Sz2 m' n') (Image arr) =
+  -- Image (A.makeArray (A.getComp arr) sz' getNewPx)
   -- where
-  --   (center@(u :. _), neighborhood) = interpolationBox method
-  --   darr =
-  --     A.makeArray
-  --       (A.getComp arr)
-  --       sz'
-  --       (getNewPx (A.handleBorderIndex border sz (A.index' arr)))
-  --   warr =
-  --     A.insertWindow darr $
-  --     A.Window
-  --       { A.windowStart = center
-  --       , A.windowSize = sz - neighborhood + Sz center
-  --       , A.windowIndex = getNewPx (A.unsafeIndex arr)
-  --       , A.windowUnrollIx2 = Just u
-  --       }
   --   sz@(Sz2 m n) = A.size arr
   --   !fM = fromIntegral m' / fromIntegral m
   --   !fN = fromIntegral n' / fromIntegral n
-  --   getNewPx getOldPx (i :. j) =
+  --   getNewPx (i :. j) =
   --     interpolate
   --       method
-  --       getOldPx
+  --       (A.handleBorderIndex border sz (A.index' arr))
   --       ( (fromIntegral i + 0.5) / fM - 0.5
-  --       , (fromIntegral j + 0.5) / fN - (0.5 :: Double))
+  --       , (fromIntegral j + 0.5) / fN - 0.5)
   --   {-# INLINE getNewPx #-}
+  Image $ A.compute warr
+  where
+    (center@(u :. _), neighborhood) = interpolationBox method
+    !darr = A.makeArray (A.getComp arr) sz' (getNewPx (A.handleBorderIndex border sz (A.index' arr)))
+    !warr =
+      A.insertWindow
+        darr
+        A.Window
+          { A.windowStart = center
+          , A.windowSize = sz - neighborhood + Sz center
+          , A.windowIndex = getNewPx (A.unsafeIndex arr)
+          , A.windowUnrollIx2 = Just u
+          }
+    sz@(Sz2 m n) = A.size arr
+    !fM = fromIntegral m' / fromIntegral m
+    !fN = fromIntegral n' / fromIntegral n
+    getNewPx getOldPx (i :. j) =
+      interpolate
+        method
+        getOldPx
+        ((fromIntegral i + 0.5) / fM - 0.5, (fromIntegral j + 0.5) / fN - 0.5)
+    {-# INLINE getNewPx #-}
 {-# INLINE resize #-}
 
+-- Note: Reducing the size seems to be better performance wise with windowed array, while
+-- increasing not necesserally
 
 -- | Scale an image. Same as resize, except scaling factors are supplied
 -- instead of new dimensions.
 --
 -- @ scale 'Bilinear' 'Edge' (0.5, 2) frog == resize 'Bilinear' 'Edge' (100, 640) frog @
 --
-scale :: (ColorModel cs e, Interpolation method) =>
+scale :: (RealFloat e, ColorModel cs e, Interpolation method) =>
          method -- ^ Interpolation method to be used during scaling.
       -> Border (Pixel cs e) -- ^ Border handling strategy
       -> (Double, Double) -- ^ Positive scaling factors.
@@ -409,16 +408,15 @@ scale method border (fM, fN) img =
 ----------------------
 
 -- | Put an angle into @[0, 2*pi)@ range.
-angle0to2pi :: Double -> Double
+angle0to2pi :: RealFloat e => e -> e
 angle0to2pi !f = f - 2 * pi * floor' (f / (2 * pi))
- where  floor' :: Double -> Double
-        floor' !x = fromIntegral (floor x :: Int)
+ where  floor' !x = fromIntegral (floor x :: Int)
         {-# INLINE floor' #-}
 {-# INLINE angle0to2pi #-}
 
 
 -- | Make sure @sin' pi == 0@ instead of @sin pi == 1.2246467991473532e-16@
-sin' :: Double -> Double
+sin' :: RealFloat e => e -> e
 sin' a = if abs sinA <= _0 then 0 else sinA
   where !_0 = 10 * sin pi
         !sinA = sin a
@@ -427,7 +425,7 @@ sin' a = if abs sinA <= _0 then 0 else sinA
 
 -- | Make sure @cos' (pi/2) == 0@ instead of @cos (pi/2) == 6.123233995736766e-17@
 -- and @cos' (3*pi/2) == 0@ instead of @cos (3*pi/2) == -1.8369701987210297e-16@
-cos' :: Double -> Double
+cos' :: RealFloat e => e -> e
 cos' a = sin' (a + pi/2)
 {-# INLINE cos' #-}
 

hip/src/Graphics/Image/Processing/Interpolation.hs

@@ -25,12 +25,12 @@ class Interpolation method where
   interpolationBox :: method -> (Ix2, Sz2)
 
   -- | Construct a new pixel by using information from neighboring pixels.
-  interpolate :: (Elevator a, RealFloat a, ColorModel cs e) =>
+  interpolate :: (RealFloat e, ColorModel cs e) =>
                  method -- ^ Interpolation method
               -> (Ix2 -> Pixel cs e)
                  -- ^ Lookup function that returns a pixel at @i@th and @j@th
                  -- location.
-              -> (a, a) -- ^ Real values of @i@ and @j@ index
+              -> (e, e) -- ^ Real values of @i@ and @j@ index
               -> Pixel cs e
 
 
@@ -65,14 +65,15 @@ instance Interpolation Bilinear where
       !j0 = floor j
       !i1 = i0 + 1
       !j1 = j0 + 1
-      !iWeight = fromRealFloat (i - fromIntegral i0)
-      !jWeight = fromRealFloat (j - fromIntegral j0)
+      !iWeight = i - fromIntegral i0
+      !jWeight = j - fromIntegral j0
       !f00 = getPx (i0 :. j0)
       !f10 = getPx (i1 :. j0)
       !f01 = getPx (i0 :. j1)
       !f11 = getPx (i1 :. j1)
       !fi0 = f00 + fmap (iWeight *) (f10 - f00)
       !fi1 = f01 + fmap (iWeight *) (f11 - f01)
+  {-# INLINE interpolate #-}
 
 
 
@@ -81,15 +82,11 @@ instance Interpolation Bicubic where
   interpolationBox _ = (1 :. 1, Sz (4 :. 4))
 
   interpolate (Bicubic a) getPx (i, j) =
-      (// fromRealFloat w) <$> ( f00 + f10 + f20 + f30
-                               + f01 + f11 + f21 + f31
-                               + f02 + f12 + f22 + f32
-                               + f03 + f13 + f23 + f33 )
+      (/ w) <$> ( f00 + f10 + f20 + f30
+                + f01 + f11 + f21 + f31
+                + f02 + f12 + f22 + f32
+                + f03 + f13 + f23 + f33 )
     where
-      distX x = fromIntegral x - i
-      {-# INLINE distX #-}
-      distY y = fromIntegral y - j
-      {-# INLINE distY #-}
       a' = fromDouble a
       weight x
           | x' <= 1 = ((a' + 2) * x' - (a' + 3)) * x2' + 1
@@ -107,14 +104,14 @@ instance Interpolation Bicubic where
       !i3 = i1 + 2
       !j3 = j1 + 2
 
-      !weightX0 = weight (distX i0)
-      !weightX1 = weight (distX i1)
-      !weightX2 = weight (distX i2)
-      !weightX3 = weight (distX i3)
-      !weightY0 = weight (distY j0)
-      !weightY1 = weight (distY j1)
-      !weightY2 = weight (distY j2)
-      !weightY3 = weight (distY j3)
+      !weightX0 = weight (fromIntegral i0 - i)
+      !weightX1 = weight (fromIntegral i1 - i)
+      !weightX2 = weight (fromIntegral i2 - i)
+      !weightX3 = weight (fromIntegral i3 - i)
+      !weightY0 = weight (fromIntegral j0 - j)
+      !weightY1 = weight (fromIntegral j1 - j)
+      !weightY2 = weight (fromIntegral j2 - j)
+      !weightY3 = weight (fromIntegral j3 - j)
 
       !weightX0Y0 = weightX0 * weightY0
       !weightX1Y0 = weightX1 * weightY0
@@ -136,25 +133,25 @@ instance Interpolation Bicubic where
       !weightX2Y3 = weightX2 * weightY3
       !weightX3Y3 = weightX3 * weightY3
 
-      !f00 = (fromRealFloat weightX0Y0 *) <$> getPx (i0 :. j0)
-      !f10 = (fromRealFloat weightX1Y0 *) <$> getPx (i1 :. j0)
-      !f20 = (fromRealFloat weightX2Y0 *) <$> getPx (i2 :. j0)
-      !f30 = (fromRealFloat weightX3Y0 *) <$> getPx (i3 :. j0)
-
-      !f01 = (fromRealFloat weightX0Y1 *) <$> getPx (i0 :. j1)
-      !f11 = (fromRealFloat weightX1Y1 *) <$> getPx (i1 :. j1)
-      !f21 = (fromRealFloat weightX2Y1 *) <$> getPx (i2 :. j1)
-      !f31 = (fromRealFloat weightX3Y1 *) <$> getPx (i3 :. j1)
-
-      !f02 = (fromRealFloat weightX0Y2 *) <$> getPx (i0 :. j2)
-      !f12 = (fromRealFloat weightX1Y2 *) <$> getPx (i1 :. j2)
-      !f22 = (fromRealFloat weightX2Y2 *) <$> getPx (i2 :. j2)
-      !f32 = (fromRealFloat weightX3Y2 *) <$> getPx (i3 :. j2)
-
-      !f03 = (fromRealFloat weightX0Y3 *) <$> getPx (i0 :. j3)
-      !f13 = (fromRealFloat weightX1Y3 *) <$> getPx (i1 :. j3)
-      !f23 = (fromRealFloat weightX2Y3 *) <$> getPx (i2 :. j3)
-      !f33 = (fromRealFloat weightX3Y3 *) <$> getPx (i3 :. j3)
+      !f00 = (weightX0Y0 *) <$> getPx (i0 :. j0)
+      !f10 = (weightX1Y0 *) <$> getPx (i1 :. j0)
+      !f20 = (weightX2Y0 *) <$> getPx (i2 :. j0)
+      !f30 = (weightX3Y0 *) <$> getPx (i3 :. j0)
+
+      !f01 = (weightX0Y1 *) <$> getPx (i0 :. j1)
+      !f11 = (weightX1Y1 *) <$> getPx (i1 :. j1)
+      !f21 = (weightX2Y1 *) <$> getPx (i2 :. j1)
+      !f31 = (weightX3Y1 *) <$> getPx (i3 :. j1)
+
+      !f02 = (weightX0Y2 *) <$> getPx (i0 :. j2)
+      !f12 = (weightX1Y2 *) <$> getPx (i1 :. j2)
+      !f22 = (weightX2Y2 *) <$> getPx (i2 :. j2)
+      !f32 = (weightX3Y2 *) <$> getPx (i3 :. j2)
+
+      !f03 = (weightX0Y3 *) <$> getPx (i0 :. j3)
+      !f13 = (weightX1Y3 *) <$> getPx (i1 :. j3)
+      !f23 = (weightX2Y3 *) <$> getPx (i2 :. j3)
+      !f33 = (weightX3Y3 *) <$> getPx (i3 :. j3)
 
       !w = weightX0Y0 + weightX1Y0 + weightX2Y0 + weightX3Y0
          + weightX0Y1 + weightX1Y1 + weightX2Y1 + weightX3Y1