Bicubic intepolation (#37)

Bicubic interpolation: * Inline interpolation * Choose parameter from -0.5 to - 1 * Optimize calculation
lehins · May 3, 2020 · 9debfc3 · 9debfc3
1 parent ee1a68b
commit 9debfc3
Show file tree

Hide file tree

Showing 2 changed files with 96 additions and 2 deletions.
diff --git a/src/Graphics/Image/Processing/Interpolation.hs b/src/Graphics/Image/Processing/Interpolation.hs
@@ -10,7 +10,7 @@
 -- Portability : non-portable
 --
 module Graphics.Image.Processing.Interpolation (
-  Interpolation(..), Nearest(..), Bilinear(..)
+  Interpolation(..), Nearest(..), Bilinear(..), Bicubic(..)
   ) where
 
 import Graphics.Image.Interface
@@ -38,6 +38,11 @@ data Nearest = Nearest deriving Show
 data Bilinear = Bilinear deriving Show
 
 
+-- | Bicubic interpolation method.
+--   The parameter is usually set from -0.5 to -1.0.
+data Bicubic = Bicubic Double deriving Show
+
+
 instance Interpolation Nearest where
 
   interpolate Nearest border !sz getPx !(i, j) =
@@ -61,3 +66,85 @@ instance Interpolation Bilinear where
     !fi0 = f00 + iPx*(f10-f00)
     !fi1 = f01 + iPx*(f11-f01)
   {-# INLINE interpolate #-}
+
+
+instance Interpolation Bicubic where
+
+  interpolate (Bicubic a) border !sz getPx !(i, j) = ( f00 + f10 + f20 + f30
+                                                     + f01 + f11 + f21 + f31
+                                                     + f02 + f12 + f22 + f32
+                                                     + f03 + f13 + f23 + f33
+                                                     ) * w  where
+    getPx' = handleBorderIndex border sz getPx
+    {-# INLINE getPx' #-}
+    distX x = fromIntegral x - i
+    {-# INLINE distX #-}
+    distY y = fromIntegral y - j
+    {-# INLINE distY #-}
+    weight x
+        | x' <= 1 = (a + 2) * x' ** 3 - (a + 3) * x' ** 2 + 1
+        | x' <  2 = a * x' ** 3 - 5 * a * x' ** 2 + 8 * a * x' - 4 * a
+        | otherwise = 0
+        where x' = abs x
+    {-# INLINE weight #-}
+    !(i1, j1) = (floor i, floor j)
+    !(i0, j0) = (i1 - 1, j1 - 1)
+    !(i2, j2) = (i1 + 1, j1 + 1)
+    !(i3, j3) = (i1 + 2, j1 + 2)
+
+    !weightX0 = weight (distX i0)
+    !weightY0 = weight (distY i0)
+    !weightX1 = weight (distX i1)
+    !weightY1 = weight (distY i1)
+    !weightX2 = weight (distX i2)
+    !weightY2 = weight (distY i2)
+    !weightX3 = weight (distX i3)
+    !weightY3 = weight (distY i3)
+
+    !weightX0Y0 = weightX0 * weightY0
+    !weightX1Y0 = weightX1 * weightY0
+    !weightX2Y0 = weightX2 * weightY0
+    !weightX3Y0 = weightX3 * weightY0
+
+    !weightX0Y1 = weightX0 * weightY1
+    !weightX1Y1 = weightX1 * weightY1
+    !weightX2Y1 = weightX2 * weightY1
+    !weightX3Y1 = weightX3 * weightY1
+
+    !weightX0Y2 = weightX0 * weightY2
+    !weightX1Y2 = weightX1 * weightY2
+    !weightX2Y2 = weightX2 * weightY2
+    !weightX3Y2 = weightX3 * weightY2
+
+    !weightX0Y3 = weightX0 * weightY3
+    !weightX1Y3 = weightX1 * weightY3
+    !weightX2Y3 = weightX2 * weightY3
+    !weightX3Y3 = weightX3 * weightY3
+
+    !f00 = getPx' (i0, j0) * promote (fromDouble weightX0Y0)
+    !f10 = getPx' (i1, j0) * promote (fromDouble weightX1Y0)
+    !f20 = getPx' (i2, j0) * promote (fromDouble weightX2Y0)
+    !f30 = getPx' (i3, j0) * promote (fromDouble weightX3Y0)
+
+    !f01 = getPx' (i0, j1) * promote (fromDouble weightX0Y1)
+    !f11 = getPx' (i1, j1) * promote (fromDouble weightX1Y1)
+    !f21 = getPx' (i2, j1) * promote (fromDouble weightX2Y1)
+    !f31 = getPx' (i3, j1) * promote (fromDouble weightX3Y1)
+
+    !f02 = getPx' (i0, j2) * promote (fromDouble weightX0Y2)
+    !f12 = getPx' (i1, j2) * promote (fromDouble weightX1Y2)
+    !f22 = getPx' (i2, j2) * promote (fromDouble weightX2Y2)
+    !f32 = getPx' (i3, j2) * promote (fromDouble weightX3Y2)
+
+    !f03 = getPx' (i0, j3) * promote (fromDouble weightX0Y3)
+    !f13 = getPx' (i1, j3) * promote (fromDouble weightX1Y3)
+    !f23 = getPx' (i2, j3) * promote (fromDouble weightX2Y3)
+    !f33 = getPx' (i3, j3) * promote (fromDouble weightX3Y3)
+
+    !w = promote . fromDouble . (1 /) $
+          weightX0Y0 + weightX1Y0 + weightX2Y0 + weightX3Y0
+        + weightX0Y1 + weightX1Y1 + weightX2Y1 + weightX3Y1
+        + weightX0Y2 + weightX1Y2 + weightX2Y2 + weightX3Y2
+        + weightX0Y3 + weightX1Y3 + weightX2Y3 + weightX3Y3
+  {-# INLINE interpolate #-}
+
diff --git a/tests/Graphics/Image/ProcessingSpec.hs b/tests/Graphics/Image/ProcessingSpec.hs
@@ -13,17 +13,20 @@ import Graphics.Image.InterfaceSpec (translateWrap, dummyImage10x20)
 data Interpol
   = I1 Nearest
   | I2 Bilinear
+  | I3 Bicubic
 
 instance Show Interpol where
   show (I1 i) = "I1 " ++ show i
   show (I2 i) = "I2 " ++ show i
+  show (I3 i) = "I3 " ++ show i
 
 instance Arbitrary Interpol where
   arbitrary = do
     ix <- arbitrary
-    case ix `mod` (2 :: Int) of
+    case ix `mod` (3 :: Int) of
       0 -> return $ I1 Nearest
       1 -> return $ I2 Bilinear
+      2 -> I3 . Bicubic <$> choose (-1, -0.5)
       _ -> error $ "Unknown interpolation: " ++ show ix
 
 
@@ -74,18 +77,22 @@ prop_concatRotate img =
 prop_rotate90 :: Interpol -> Border (Pixel RGB Double) -> Image VU RGB Double -> Bool
 prop_rotate90 (I1 i) border img = rotate90 img == rotate i border (pi/2) img
 prop_rotate90 (I2 i) border img = rotate90 img == rotate i border (pi/2) img
+prop_rotate90 (I3 i) border img = rotate90 img == rotate i border (pi/2) img
 
 prop_rotate180 :: Interpol -> Border (Pixel RGB Double) -> Image VU RGB Double -> Bool
 prop_rotate180 (I1 i) border img = rotate180 img == rotate i border pi img
 prop_rotate180 (I2 i) border img = rotate180 img == rotate i border pi img
+prop_rotate180 (I3 i) border img = rotate180 img == rotate i border pi img
 
 prop_rotate270 :: Interpol -> Border (Pixel RGB Double) -> Image VU RGB Double -> Bool
 prop_rotate270 (I1 i) border img = rotate270 img == rotate i border (3*pi/2) img
 prop_rotate270 (I2 i) border img = rotate270 img == rotate i border (3*pi/2) img
+prop_rotate270 (I3 i) border img = rotate270 img == rotate i border (3*pi/2) img
 
 prop_rotate360 :: Interpol -> Border (Pixel RGB Double) -> Image VU RGB Double -> Bool
 prop_rotate360 (I1 i) border img = (rotate270 . rotate90) img == rotate i border (2*pi) img
 prop_rotate360 (I2 i) border img = (rotate270 . rotate90) img == rotate i border (2*pi) img
+prop_rotate360 (I3 i) border img = (rotate270 . rotate90) img == rotate i border (2*pi) img
 
 
 struct :: Image VS X Bit