Transparency aware photo image rotation

EG - Dealing with image transformations, I stumbled upon BLT's Picture image type, which has two interesting features: rotation and resampling (aka scaling). Here is a translation to pure Tcl of its image rotation code, which deals with the alpha channel and performs a decent antialiasing as described in Rotating a Tk Photo Image. The code in turn is derived from the Leptonica image library.

Note that this code doesn't handle the special case of 90°/180°/270° degree rotations, which can be done way faster using other methods.

package require Tk

namespace eval imgrot {
    namespace import ::tcl::mathop::+
}

proc imgrot::Getrotatedsize {w h angle} {
    set angle [expr {$angle * acos(-1) / 180.0}]
    set sina  [expr {sin($angle)}]
    set cosa  [expr {cos($angle)}]

    # Set the four corners of the rectangle whose center is the origin.
    #  0---1
    #  |   |
    #  3---2
    #

    set corner(1.x) [set corner(2.x) [expr {$w * 0.5}]]
    set corner(0.x) [set corner(3.x) [expr {-$corner(1.x)}]];
    set corner(2.y) [set corner(3.y) [expr {$h * 0.5}]]
    set corner(0.y) [set corner(1.y) [expr {-$corner(2.y)}]]

    set xMax [set yMax 0.0]]

    # Rotate the four corners and find the maximum X and Y coordinates
    for {set i 0} {$i < 4} {incr i} {

        set x [expr {($corner($i.x) * $cosa) - ($corner($i.y) * $sina)}]
        set y [expr {($corner($i.x) * $sina) + ($corner($i.y) * $cosa)}]
        if {$x > $xMax} {
            set xMax $x
        }
        if {$y > $yMax} {
            set yMax $y
        }
    }
    # By symmetry, the width and height of the bounding box are twice the
    # maximum x and y coordinates.
    return [list [expr {int(floor(2 * $xMax + 0.5))}] \
        [expr {int(floor(2 * $yMax + 0.5))}]]
}

proc imgrot::UCLAMP {v} {
    expr {$v > 255 ? 255 : ($v < 0 ? 0 : $v)}
}

# use the bg parameter to fill the areas not covered by the source image
proc imgrot::rotatebyareamapping {srcimg dstimg angle {bg "#ffffff00"}} {
    set srcw [image width  $srcimg]
    set srch [image height $srcimg]
    lassign [Getrotatedsize $srcw $srch $angle] rotWidth rotHeight
    $dstimg configure -width $rotWidth -height $rotHeight
    $dstimg blank

    set srcCx [expr {$srcw / 2}]
    set srcCy [expr {$srch / 2}]
    set destCx [expr {$rotWidth / 2}]
    set destCy [expr {$rotHeight / 2}]
    set wm2 [expr {$srcw - 2}]
    set hm2 [expr {$srch - 2}]

    set radians [expr {-$angle * acos(-1) / 180.0}]
    set sinTheta [expr {16.0 * sin($radians)}]
    set cosTheta [expr {16.0 * cos($radians)}]

    for {set y  0} {$y < $rotHeight} {incr y} {
        set deltaY [expr {$destCy - $y}]

        for {set x 0} {$x < $rotWidth} {incr x} {
            set deltaX [expr {$destCx - $x}]
            set xpm [expr {int((-$deltaX * $cosTheta) - ($deltaY * $sinTheta))}]
            set ypm [expr {int((-$deltaY * $cosTheta) + ($deltaX * $sinTheta))}]
            set srcX [expr {$srcCx + ($xpm >> 4)}]
            set srcY [expr {$srcCy + ($ypm >> 4)}]
            set xf [expr {$xpm & 0x0f}]
            set yf [expr {$ypm & 0x0f}]

            # If outside of the source image, use the default color
            if {($srcX < 0) || ($srcY < 0) || ($srcX > $wm2) || ($srcY > $hm2)} {
                $dstimg put $bg -to $x $y
                continue
            }
            # do area weighting.  Without this, we would
            # simply do:
            #  *(lined + x) = *(lines + srcX);
            # which is faster but gives lousy results!
            #
            lassign [$srcimg get $srcX $srcY -withalpha] p00r p00g p00b p00a
            lassign [$srcimg get [+ $srcX 1] $srcY -withalpha] p01r p01g p01b p01a
            lassign [$srcimg get $srcX [+ $srcY 1] -withalpha] p10r p10g p10b p10a
            lassign [$srcimg get [+ $srcX 1] [+ $srcY 1] -withalpha] p11r p11g p11b p11a

            set r [expr {
                 ((16 - $xf) * (16 - $yf) * $p00r +
                 $xf * (16 - $yf) * $p01r + (16 - $xf) * $yf * $p10r +
                 $xf * $yf * $p11r + 128) / 256
            }]
            set g [expr {
                 ((16 - $xf) * (16 - $yf) * $p00g +
                 $xf * (16 - $yf) * $p01g + (16 - $xf) * $yf * $p10g +
                 $xf * $yf * $p11g + 128) / 256
            }]
            set b [expr {
                 ((16 - $xf) * (16 - $yf) * $p00b +
                 $xf * (16 - $yf) * $p01b + (16 - $xf) * $yf * $p10b +
                 $xf * $yf * $p11b + 128) / 256
            }]
            set a [expr {
                ((16 - $xf) * (16 - $yf) * $p00a +
                $xf * (16 - $yf) * $p01a + (16 - $xf) * $yf * $p10a +
                $xf * $yf * $p11a + 128) / 256
            }]

            $dstimg put [format {#%02X%02X%02X%02X} \
                    [UCLAMP $r] [UCLAMP $g] [UCLAMP $b] [UCLAMP $a]] \
                -to $x $y
        }
    }
}

# Demo code
# use the feather image from the library/demos/images Tk source directory
set srcimg [image create photo \
        -file [file join $env(HOME) src tk library demos images Tk_feather.png]]
set dstimg [image create photo]
set angle 0.0
label .langle -text Angle
spinbox .angle -textvariable angle -from 0.0 -to 359 -increment 1 \
    -validatecommand {string is double -strict %P} \
    -validate all
button .b -text Rotate -width 15 -command {
    imgrot::rotatebyareamapping $srcimg $dstimg $angle
}
label .src -image $srcimg
label .dst -image $dstimg

grid .langle .angle .b -padx 3 -pady 3
grid .src - - -pady 3
grid .dst - - -pady 3