Now, it looks like you're doing it right. I expected this, see
https://github.com/timholy/Grid.jl/pull/38. This is part of what I meant by
"refactoring" :). However, for image interpolation even further savings beyond
that pull request are possible: for example, you only need one call to floor
per pixel, because you can nest the index operations.
Here's a reasonably well-optimized prototype (a flexible implementation will
take longer to write, but should not cost performance). Note there may be some
additional optimizations possible.
First, the results:
julia> include("/tmp/resize.jl");
Cairo:
elapsed time: 0.098905264 seconds (53896192 bytes allocated, 38.25% gc time)
Julia:
elapsed time: 0.034537582 seconds (6340816 bytes allocated)
Julia is ~3 times faster.
Here's the code (I just copied your Cairo implementation):
using Images, Cairo
function imresize_julia!(resized, original)
scale1 = (size(original,1)-1)/(size(resized,1)-0.999f0)
scale2 = (size(original,2)-1)/(size(resized,2)-0.999f0)
for jr = 0:size(resized,2)-1
jo = scale2*jr
ijo = itrunc(jo)
fjo = jo - oftype(jo, ijo)
@inbounds for ir = 0:size(resized,1)-1
io = scale1*ir
iio = itrunc(io)
fio = io - oftype(io, iio)
tmp = (1-fio)*((1-fjo)*original[iio+1,ijo+1] +
fjo*original[iio+1,ijo+2])
+ fio*((1-fjo)*original[iio+2,ijo+1] +
fjo*original[iio+2,ijo+2])
resized[ir+1,jr+1] = convertsafely(eltype(resized), tmp)
end
end
resized
end
imresize_julia(original, new_size) = imresize_julia!(similar(original,
new_size), original)
convertsafely{T<:FloatingPoint}(::Type{T}, val) = convert(T, val)
convertsafely{T<:Integer}(::Type{T}, val::Integer) = convert(T, val)
convertsafely{T<:Integer}(::Type{T}, val::FloatingPoint) = itrunc(T,
val+oftype(val, 0.5))
function imresize_cairo(dat::Array{Uint32, 2}, new_size::(Int, Int))
cs = CairoImageSurface(dat, 0)
new_dat = zeros(Uint32, new_size)
new_cs = CairoImageSurface(new_dat, 0)
pat = CairoPattern(cs)
pattern_set_filter(pat, Cairo.FILTER_BILINEAR)
c = CairoContext(new_cs)
h, w = size(dat)
new_h, new_w = new_size
scale(c, new_h / h, new_w / w)
set_source(c, pat)
paint(c)
return new_cs.data
end
img = rand(0x00:0xff, 774, 512)
new_size = (3096, 2048)
imresize_cairo(convert(Array{Uint32}, img), new_size)
println("Cairo:")
@time imresize_cairo(convert(Array{Uint32}, img), new_size)
imresize_julia(img, new_size)
println("Julia:")
@time imresize_julia(img, new_size)
More advantages of doing this in Julia rather than through Cairo:
- No binary dependencies (as you stated earlier)
- Cairo is basically 8-bit, which is inadequate for many applications. Cairo's
Uint32s are actually for encoding color, and internally they get treated like
4 Uint8s---you can't get 16-bit dynamic range, for example.
- For Cairo you'd have to convert images with different datatypes into Uint32.
Measuring the performance of Cairo should include the cost of conversion. With
Julia, we can write versions that work for any input. For example, try a
Float32 image, you'll see the julia version is even faster that what I showed
above.
- We can separately implement algorithms for grayscale and color. Part of the
reason this is three times faster than Cairo's is that Cairo is basically
doing three times the work.
--Tim
On Friday, August 08, 2014 12:54:24 AM Andrei wrote:
> I've implemented Grid-based version, but it seems much slower than both -
> Cairo and mine. Am I using Grid the wrong way? Or does Grid just take care
> of more use cases?
>
> -----
>
> function imresize_grid{T <: Float64}(dat::Array{T, 2}, new_size::(Int, Int))
> new_dat = similar(dat, new_size)
> h, w = size(dat)
> new_h, new_w = new_size
> grid = InterpGrid(dat, BCnil, InterpLinear)
> for new_j=1:new_w, new_i=1:new_h
> # coordinates in original image
> i = new_i * h / new_h
> j = new_j * w / new_w
> new_dat[new_i, new_j] = grid[i, j]
> end
> return new_dat
> end
>
> ----
>
> And here are performance tests (imf64 is an Array{Float64, 2}, while imu32
> is the same image converted to Array{Uint32, 2}):
>
>
> julia> @time for i=1:1000 imresize(imf64, (128, 128)) end
> elapsed time: 0.697520586 seconds (131184000 bytes allocated, 15.41% gc
> time)
>
> julia> @time for i=1:1000 imresize_cairo(imu32, (128, 128)) end
> elapsed time: 0.288527903 seconds (393616000 bytes allocated, 24.78% gc
> time)
>
> julia> @time for i=1:1000 imresize_grid(imf64, (128, 128)) end
> elapsed time: 4.759103481 seconds (1443800000 bytes allocated, 14.29% gc
> time)
>
> On Thu, Aug 7, 2014 at 4:21 AM, Tim Holy <[email protected]> wrote:
> > I bet we can meet or beat Cairo; I've never found Cairo to be terribly
> > fast.
> > But of course one doesn't know until one tries.
> >
> > Given that the main issue is really interpolation, it (and some of the
> > refactoring that's currently under discussion) seems relevant. There's
> > some
> > work brewing on a new version of Grid, and being able to do this
> > performantly
> > and easily can become a goal there. Grid already has more interpolation
> > options than Cairo does.
> >
> > --Tim
> >
> > On Thursday, August 07, 2014 01:41:39 AM Andrei wrote:
> > > @Andreas: thanks and sorry for late reply.
> > >
> > > Using your code, I created testable function that works almost 10x times
> > > faster then my naive implementation! Right now, though, it rotates image
> >
> > 90
> >
> > > degrees to the left for some reason I don't really understand, but I
> > > believe it's fixable.
> > >
> > > Here's a code for (not fully correct) function:
> > >
> > > function imresize_cairo(dat::Array{Uint32, 2}, new_size::(Int, Int))
> > >
> > > cs = CairoImageSurface(dat, 0)
> > > new_dat = zeros(Uint32, new_size)
> > > new_cs = CairoImageSurface(new_dat, 0)
> > > pat = CairoPattern(cs)
> > > pattern_set_filter(pat, Cairo.FILTER_BILINEAR)
> > > c = CairoContext(new_cs)
> > > h, w = size(dat)
> > > new_h, new_w = new_size
> > > scale(c, new_h / h, new_w / w)
> > > set_source(c, pat)
> > > paint(c)
> > > return new_cs.data
> > >
> > > end
> > >
> > >
> > > I'm pretty satisfied with this version (or what it should turn into)
> >
> > both -
> >
> > > because of its speed and possibility to use different interpolation
> > > schemas. Though, one disadvantage of using it in public package like
> > > Images.jl is that it brings additional dependency, and we all know how
> > > annoying these dependencies may be sometimes. We could implement 2
> >
> > versions
> >
> > > - one with Cairo and fast and another simple and slow - and then load
> > > appropriate function via macros. But since simple version doesn't
> > > support
> > > interpolation other than bilinear, we will have to restrict Cairo
> > > version
> > > to this type too to keep same interface. Or we can just 2 different
> > > functions, but it may be pretty confusing for new users.
> > >
> > > Some opinions on a better approach would be helpful here.
> > >
> > >
> > >
> > >
> > > On Tue, Aug 5, 2014 at 11:24 AM, Andreas Lobinger <[email protected]>
> > >
> > > wrote:
> > > > An example looks like this:
> > > >
> > > > using Cairo
> > > >
> > > > # prepapration, an image of 2x2 pixels, an target of 8x8 pixels
> > > >
> > > > d = [8 4; 2 1];
> > > > d0 = reinterpret(Uint32,d);
> > > > cs0 = Cairo.CairoImageSurface(d0,0); # 0 as FORMAT_ARGB32
> > > > show(cs0.data)
> > > > print('\n')
> > > >
> > > > d2 = zeros(Uint32,8,8);
> > > > cs2 = Cairo.CairoImageSurface(d2,0);
> > > >
> > > > p = Cairo.CairoPattern(cs0);
> > > >
> > > > # just scale by using FILTER_NEAREST
> > > >
> > > > pattern_set_filter(p,Cairo.FILTER_NEAREST);
> > > > c = Cairo.CairoContext(cs2);
> > > > scale(c,4,4);
> > > > set_source(c,p);
> > > > paint(c);
> > > > show(cs2.data)
> > > > print('\n')
> > > >
> > > > # now interpolate
> > > >
> > > > d3 = zeros(Uint32,8,8);
> > > > cs3 = Cairo.CairoImageSurface(d3,0);
> > > >
> > > > pattern_set_filter(p,Cairo.FILTER_GOOD);
> > > > c = Cairo.CairoContext(cs3);
> > > > scale(c,4,4);
> > > > set_source(c,p);
> > > > paint(c);
> > > > show(cs3.data)
> > > > print('\n')
> > > >
> > > > and gives on my command line:
> > > >
> > > > julia> include("ca.jl")
> > > > Uint32[8 2
> > > >
> > > > 4 1]
> > > >
> > > > Uint32[8 8 8 8 2 2 2 2
> > > >
> > > > 8 8 8 8 2 2 2 2
> > > > 8 8 8 8 2 2 2 2
> > > > 8 8 8 8 2 2 2 2
> > > > 4 4 4 4 1 1 1 1
> > > > 4 4 4 4 1 1 1 1
> > > > 4 4 4 4 1 1 1 1
> > > > 4 4 4 4 1 1 1 1]
> > > >
> > > > Uint32[3 4 4 3 2 1 1 0
> > > >
> > > > 4 6 6 5 3 2 1 1
> > > > 4 6 6 5 3 2 1 1
> > > > 4 5 5 4 3 2 1 1
> > > > 3 4 4 3 2 1 1 0
> > > > 2 3 4 3 2 1 0 0
> > > > 2 3 3 2 1 1 0 0
> > > > 1 2 2 1 1 0 0 0]
> > > >
> > > > But this is just how to handle the data. For interpreting the pixels
> > > > correctly you need to read little bit about the colormodel and pixel
> > > > formats.
> > > >
> > > >
> > > > Wishing a happy day,
> > > >
> > > > Andreas
