def op_l0smoothing_v {
label {L0 Smoothing, HSV, V-only}
setup_image {
showsmooth [base]
}
setup {
proc showsmooth {i} {
# Do the smoothing separately on each channel. While this
# is not how the original performs the operation it is
# easier to implement. I will have to see how much this
# affects the result.
lassign [crimp split [crimp convert 2hsv $i]] h s v
set v [l0smooth $v 0.005]
show_image [crimp montage vertical \
[crimp convert 2rgb [crimp join 2hsv $h $s $v]] \
[crimp convert 2rgb $i]]
return
}
namespace eval l0smooth {
variable DX [crimp kernel fpmake {{ 0 1 -1}} 0]
variable DY [crimp kernel transpose $DX]
# DX, DY :: float, convolution kernel, spatial
variable DXI [crimp read tcl float {{1 -1}}]
variable DYI [crimp read tcl float {{1} {-1}}]
# DXI, DYI :: float, convolution kernel, spatial, FFT inputs.
variable tofloat [expr {1./255}]
variable toint 255.0
}
proc l0smooth::Setup {} {
variable DX
variable DY
variable DXI
variable DYI
variable tofloat
upvar 1 lambda lambda ; # :: scalar float
upvar 1 beta0 beta0 ; # :: scalar float
upvar 1 image image ; # :: img grey8 -> float, spatial
upvar 1 result result ; # :: img float, spatial
upvar 1 beta beta ; # :: scalar float
upvar 1 count count ; # :: scalar int
upvar 1 f1 f1 ; # :: img complex, frequency
upvar 1 fi fi ; # :: img complex, frequency
upvar 1 cfdx cfdx ; # :: img complex, frequency
upvar 1 cfdy cfdy ; # :: img complex, frequency
upvar 1 denom denom ; # :: img float, frequency
# beta0 default = 2 lambda
if {$beta0 eq {}} {
set beta0 [expr {2* $lambda}]
}
set image [crimp::scale_float \
[crimp convert 2float $image] $tofloat]
# image :: img float, range [0...1]
# loop initialization
set result $image
set beta $beta0
set count 0
lassign [crimp dimensions $image] w h
# Constant parts of the loop moved upfront.
set f1 [crimp convert 2complex [crimp blank float $w $h 1.0]]
# f1 :: img complex, frequency
# This is the FFT of the delta function. That is, a single
# 1 in the top left corner. Interestingly, the result is
# an image with an 1 in all pixels. So we spare us the FFT
# and generate this directly.
set fi [crimp fft forward [crimp convert 2complex $image]]
# fi :: img complex, frequency
# FFT of the convolution kernels ...
# It is easier to construct them directly than to pull the
# images out of the kernels, and then shift them properly
# during the expansion. (The center 1 must be placed in
# the top left corner).
# We assume that DX/DY are smaller than the image in all
# respects. We will be in trouble should the input image
# get smaller than 3x3.
log dxi=[crimp::dimensions $DXI]
log dyi=[crimp::dimensions $DYI]
set dx [crimp expand const $DXI 0 0 [expr {$w-2}] [expr {$h-1}] 0.0]
set dy [crimp expand const $DYI 0 0 [expr {$w-1}] [expr {$h-2}] 0.0]
# dx, dy :: img float, spatial
log dx=[crimp::dimensions $dx]
log dy=[crimp::dimensions $dy]
set fdx [crimp fft forward [crimp convert 2complex $dx]]
set fdy [crimp fft forward [crimp convert 2complex $dy]]
# fdx, fdy :: img complex, frequency
# Conjugates. See Solve8.
set cfdx [crimp complex conjugate $fdx]
set cfdy [crimp complex conjugate $fdy]
# cfdx, cfdy :: img complex, frequency
# The main part of the denominator used in Solve8 uses
# multiplication of numbers with their own complex
# conjugate (fdX, cfdX). That comes down to the squared
# magnitude, which we can compute easier. We pre-convert
# the result to complex as that allows us to avoid a few
# casts (conversions) later in the loop.
set denom [crimp convert 2complex \
[crimp add \
[crimp::sqmagnitude_fpcomplex $fdx] \
[crimp::sqmagnitude_fpcomplex $fdy]]]
# denom :: img complex, frequency
return
}
proc l0smooth::Solve12 {} {
variable DX
variable DY
upvar 1 lambda lambda ; # :: scalar float
upvar 1 beta beta ; # :: scalar float
upvar 1 h h ; # :: img float, spatial
upvar 1 v v ; # :: img float, spatial
upvar 1 result result ; # :: img float, spatial
log \tEq/12;update
# Equation 12: Solve for (new) h, v using the last round
# result.
# TODO - The paper defines gradient as difference between
# adjacent pixels => Convolve [1,-1], and transposed.
# Question is, can we do the fast gaussian 01/10 commands
# here, with a fixed sigma ?
# TODO: Gradient for RGB it seems to be euclidean distance
# between pixel values treated as 3-vectors. Investigate
# what happens if we solve this either separately in each
# channel, or just in luma instead of using this RGB
# distance.
# result :: img float, spatial
set dx [crimp filter convolve $result $DX]
set dy [crimp filter convolve $result $DY]
# dx, dy :: img float, spatial
set threshold [crimp add [crimp square $dx] [crimp square $dy]]
# threshold :: img float, spatial
# h = << 0, for threshold <= lambda/beta, and dx otherwise >>
# v = << 0, for threshold <= lambda/beta, and dy otherwise >>
# mult = << 0, for threshold <= lambda/beta, and 1 otherwise >>
# = (T <= t) ? 0 : 1
# = (T > t) ? 1 : 0
# = above (T, t)
set threshold [crimp threshold global above \
$threshold \
[expr {$lambda / $beta}]]
# threshold :: img float, semi-boolean, spatial
# At last, multiplication of this with dx, dy
# generates the sought for h, v, removing the unwanted
# gradients and leaving the others untouched.
set h [crimp multiply $dx $threshold]
set v [crimp multiply $dy $threshold]
# h, v :: img float, spatial
return
}
proc l0smooth::Solve8 {} {
upvar 1 beta beta ; # :: scalar float
upvar 1 result result ; # :: img float, spatial
upvar 1 h h ; # :: img float, spatial
upvar 1 v v ; # :: img float, spatial
upvar 1 f1 f1 ; # :: img complex, frequency
upvar 1 fi fi ; # :: img complex, frequency
upvar 1 cfdx cfdx ; # :: img complex, frequency
upvar 1 cfdy cfdy ; # :: img complex, frequency
upvar 1 denom denom ; # :: img float, frequency
log \tEq/8;update
# Equation 8. Solve for a new result using the new h and v.
# NOTE: The FFTs are done in the complex domain.
# NOTE 2: Investigate if we can avoid real/complex/real
# conversions in the loop and push this before and after.
set fh [crimp fft forward [crimp convert 2complex $h]]
set fv [crimp fft forward [crimp convert 2complex $v]]
# fh, fv :: img complex, frequency
# numerator = fi + beta * (fdx* . fh + fdy* . fv)
# denominator = f1 + beta * (fdx* . fdx + fdy* . fdy)
# denom = ~~~~~~~~~~~~~~~~~~~~~~~
# + = element-wise addition
# . = element-wise multiplication
# * = element-wise complex conjugate
set numer [crimp add \
[crimp multiply $cfdx $fh] \
[crimp multiply $cfdy $fv]]
set num [crimp add $fi [crimp::scale_fpcomplex $numer $beta]]
set den [crimp add $f1 [crimp::scale_fpcomplex $denom $beta]]
# num :: img complex, frequency
# den :: img complex, frequency
set result [crimp complex real \
[crimp fft backward \
[crimp divide $num $den]]]
# result :: img float, spatial
return
}
proc l0smooth {image lambda {beta0 {}} {betamax 10000} {kappa 2}} {
variable l0smooth::toint
# beta0, betamax, kappa => defaults.
# image type -> float.
show_image $image
log "Setup l=$lambda, k=$kappa, bmax=$betamax";update
l0smooth::Setup
while {$beta < $betamax} {
log "Round $count @ $beta/$betamax"
l0smooth::Solve12
l0smooth::Solve8
# Show interim results.
show_image [crimp convert 2grey8 \
[crimp::scale_float $result $toint]]
# Prepare for next round.
set beta [expr {$beta * $kappa}]
incr count
}
log Done
return [crimp convert 2grey8 \
[crimp::scale_float $result $toint]]
}
}
}