You can't directly map the input image's pixels to the output image. If you do, some pixels don't map to integral positions so the result has holes and color discontinuities.

The key is to map pixels in the output image back to the positions they should have come from in the input image. For output position (x1, y1), you get an input position (x0, y0) where x0 and y0 are not necessarily integers. You then use a weighted average of the pixels surrounding (x1, y1) in the input image to determine the color of the output pixel.

The following code warps an image. It loops over the output image's pixels and calls subroutine MapPixel to map back to the input image position (x0, y0). It then calculates a weighted average for the result pixel's color.

' Transform the image.
Private Sub TransformImage(ByVal bm_src As Bitmap, ByVal _
    bm_dest As Bitmap, ByVal warp_op As WarpOperations)
    ' Find image information.
    Dim xmid As Double = bm_dest.Width / 2
    Dim ymid As Double = bm_dest.Height / 2
    Dim rmax As Double = bm_dest.Width * 0.75

    Dim ix_max As Integer = bm_src.Width - 2
    Dim iy_max As Integer = bm_src.Height - 2

    ' Generate a result for each output pixel.
    Dim x0 As Double
    Dim y0 As Double
    For y1 As Integer = 0 To bm_dest.Height - 1
        For x1 As Integer = 0 To bm_dest.Width - 1
            ' Map back to the source image.
            MapPixel(warp_op, xmid, ymid, rmax, x1, y1, x0, _
                y0)

            ' Interpolate to get the result pixel's value.
            ' Find the next smaller integral position.
            Dim ix0 As Integer = CInt(Int(x0))
            Dim iy0 As Integer = CInt(Int(y0))

            ' See if this is out of bounds.
            If (ix0 < 0) Or (ix0 > ix_max) Or _
               (iy0 < 0) Or (iy0 > iy_max) _
            Then
                ' The point is outside the image. Use white.
                bm_dest.SetPixel(x1, y1, Color.White)
            Else
                ' The point lies within the image.
                ' Calculate its value.
                Dim dx0 As Double = x0 - ix0
                Dim dy0 As Double = y0 - iy0
                Dim dx1 As Double = 1 - dx0
                Dim dy1 As Double = 1 - dy0

                ' Get the colors of the surrounding pixels.
                Dim color00 As Color = bm_src.GetPixel(ix0, _
                    iy0)
                Dim color01 As Color = bm_src.GetPixel(ix0, _
                    iy0 + 1)
                Dim color10 As Color = bm_src.GetPixel(ix0 _
                    + 1, iy0)
                Dim color11 As Color = bm_src.GetPixel(ix0 _
                    + 1, iy0 + 1)

                ' Compute the weighted average.
                Dim r As Integer = CInt( _
                    color00.R * dx1 * dy1 + color01.R * dx1 _
                        * dy0 + _
                    color10.R * dx0 * dy1 + color11.R * dx0 _
                        * dy0 _
                    )
                Dim g As Integer = CInt( _
                    color00.G * dx1 * dy1 + color01.G * dx1 _
                        * dy0 + _
                    color10.G * dx0 * dy1 + color11.G * dx0 _
                        * dy0 _
                    )
                Dim b As Integer = CInt( _
                    color00.B * dx1 * dy1 + color01.B * dx1 _
                        * dy0 + _
                    color10.B * dx0 * dy1 + color11.B * dx0 _
                        * dy0 _
                    )
                bm_dest.SetPixel(x1, y1, _
                    Color.FromArgb(255, r, g, b))
            End If
        Next x1
    Next y1
End Sub

' Map the output pixel (x1, y1) back to the input pixel
' (x0, y0).
Private Sub MapPixel(ByVal warp_op As WarpOperations, ByVal _
    xmid As Double, ByVal ymid As Double, ByVal rmax As _
    Double, ByVal x1 As Integer, ByVal y1 As Integer, ByRef _
    x0 As Double, ByRef y0 As Double)
    Const PI_OVER_2 As Double = PI / 2
    Const K As Double = 50
    Const OFFSET As Double = -PI_OVER_2
    Dim dx As Double
    Dim dy As Double
    Dim r1 As Double
    Dim r2 As Double

    Select Case warp_op
        Case WarpOperations.opIdentity
            x0 = x1
            y0 = y1
        Case WarpOperations.opFishEye
            dx = x1 - xmid
            dy = y1 - ymid
            r1 = Sqrt(dx * dx + dy * dy)
            If r1 = 0 Then
                x0 = xmid
                y0 = ymid
            Else
                r2 = rmax / 2 * (1 / (1 - r1 / rmax) - 1)
                x0 = dx * r2 / r1 + xmid
                y0 = dy * r2 / r1 + ymid
            End If
        Case WarpOperations.opTwist
            dx = x1 - xmid
            dy = y1 - ymid
            r1 = Sqrt(dx * dx + dy * dy)
            If r1 = 0 Then
                x0 = 0
                y0 = 0
            Else
                Dim theta As Double = Atan2(dx, dy) - r1 / _
                    K - OFFSET
                x0 = r1 * Cos(theta)
                y0 = r1 * Sin(theta)
            End If
            x0 = x0 + xmid
            y0 = y0 + ymid
        Case WarpOperations.opWave
            x0 = x1
            y0 = y1 - 10 * (Sin(x1 / 50 * PI) + 1) + 5
        Case WarpOperations.opSmallTop
            dx = xmid - x1
            dx = dx * ymid * 2 / (y1 + 1)
            x0 = xmid - dx
            y0 = y1
        Case WarpOperations.opWiggles
            dx = xmid - x1
            dx = dx * (Sin(y1 / 50 * PI) / 2 + 1.5)
            x0 = xmid - dx
            y0 = y1
        Case WarpOperations.opDoubleWave
            x0 = x1 - 10 * (Sin(y1 / 50 * PI) + 1) + 5
            y0 = y1 - 10 * (Sin(x1 / 50 * PI) + 1) + 5
        Case Else ' Flip vertically and horizontally.
            x0 = 2 * xmid - x1
            y0 = 2 * ymid - y1
    End Select
End Sub