Bicubic Interpolation

''' <summary>
  ''' Bicubically resize the src bitmap to the dest bitmap
  ''' </summary>
  ''' <param name="src"></param>
  ''' <param name="dest">Destination Bitmap, should be initialized to the requested size.</param>
  ''' <remarks>Original Source For This Function - http://www.codeproject.com/KB/recipes/aforge.aspx</remarks>    
  Sub Resize(ByVal src As Imaging.WriteableBitmap, ByVal dest As Imaging.WriteableBitmap)

      Dim startTime As DateTime = Now


      Dim srcwidth As Integer = src.PixelWidth
      Dim srcHeight As Integer = src.PixelHeight

      Dim destWidth As Integer = dest.PixelWidth
      Dim destHeight As Integer = dest.PixelHeight

      Dim xFactor As Double = srcwidth / destWidth
      Dim yFactor As Double = srcHeight / destHeight

      'Coordinates of src points and  coefficients.. WTF?
      Dim ox, oy, dx, dy, k1, k2 As Double
      Dim ox1, oy1, ox2, oy2 As Integer

      'new color values
      Dim r, g, b As Double

      Dim xmax As Integer = srcwidth - 1
      Dim ymax As Integer = srcHeight - 1

      For y As Integer = 0 To destHeight - 1
          'this is getting the original Y pixel if I'm not mistaken
          oy = y * yFactor - 0.5F
          'WHY is this necessary?
          oy1 = oy
          'WTF? It's getting the FRACTIONAL difference between
          'oy and oy1?
          dy = oy - oy1


          For x As Integer = 0 To destWidth - 1
              'same WTF as with the Y values,
              'maybe will make more sense later?
              ox = x * xFactor - 0.5F
              ox1 = ox
              dx = ox - ox1

              r = 0
              g = 0
              b = 0


              'this loop is to gather the interpolated
              'values of 2 pixels surrounding the current one
              For n As Integer = -1 To 2
                  'this gets the Y coefficient

                  k1 = BicubicKernel(dy - n)

                  'this seems to be getting the new
                  'Y pixel where the interpolated
                  'value comes from
                  oy2 = oy1 + n

                  'this is to ensure we're in the right spot
                  If oy2 < 0 Then
                      oy2 = 0
                  ElseIf oy2 > ymax Then
                      oy2 = ymax
                  End If

                  For m As Integer = -1 To 2
                      'for X coefficient
                      k2 = k1 * BicubicKernel(m - dx)

                      ox2 = ox1 + m

                      If ox2 < 0 Then
                          ox2 = 0
                      ElseIf ox2 > xmax Then
                          ox2 = xmax
                      End If

                      'get original pixel color
                      ' Dim origColor As Color = source.GetPixel(ox2, oy2)

                      'writablebitmap pixel is AARRGGBB
                      Dim srcColor As Integer = src.Pixels(ox2 + (srcwidth * oy2))


                      'set interpolated values
                      r += k2 * ((srcColor >> 16) And &HFF)
                      g += k2 * ((srcColor >> 8) And &HFF)
                      b += k2 * (srcColor And &HFF)
                  Next
              Next

              'after calculating coefficients we have our new color               
              dest.Pixels(x + (destWidth * y)) = &HFF000000 Or (r << 16) Or (g << 8) Or b


          Next
      Next
      dest.Invalidate()
      Dim endTime As DateTime = Now
      Diagnostics.Debug.WriteLine(String.Format("Time taken to squish {0}x{1} to {2}x{3} : {4} seconds", srcwidth, srcHeight, destWidth, destHeight, endTime.Subtract(startTime).TotalSeconds))

  End Sub

  Shared Function BicubicKernel(ByVal x As Double) As Double

      If x > 2.0 Then
          Return 0.0
      End If

      Dim a, b, c, d As Double
      Dim xm1 As Double = x - 1.0
      Dim xp1 As Double = x + 1.0
      Dim xp2 As Double = x + 2.0



      If xp2 <= 0 Then
          a = 0
      Else
          a = xp2 * xp2 * xp2
      End If

      If xp1 <= 0 Then
          b = 0
      Else
          b = xp1 * xp1 * xp1
      End If

      If x <= 0 Then
          c = 0
      Else
          c = x * x * x
      End If

      If xm1 <= 0 Then
          d = 0
      Else
          d = xm1 * xm1 * xm1
      End If
     
      Return (0.16666666666666666 * (a - (4.0 * b) + (6.0 * c) - (4.0 * d)))

  End Function