How to make Julia code faster than python code?

General Usage
1

I found my Juila code is way slower than python code, it keeps running for more than 2 hours and no results, so I terminated it. Python code gave results in 1 hour.

Here is the original python code:

def ImageTiler(filename, filepath, save_dir, level=-1,tile_size=224, overlap=0, limit_bounds=False, save_tiles=False, debug=False, save_empty_tiles=False, save_tile_map=True, threshold=0.95, print_tile_map=False, data_format='.mrxs'):
    """Returns all tiles of a whole slide image"""
    # file should be in .mrxs format
    slide = o.OpenSlide(filepath)
    if debug:
        print(slide, '\nfile: ', filename)
    DZG = dz.DeepZoomGenerator(slide,tile_size,overlap,limit_bounds)
    # the level is chosen
    dim = DZG.level_count + level
    # number of pixels in the current dimension(level)
    pixels = DZG.level_dimensions[dim]
    print(pixels)
    # number of tiles
    DZG_dim = DZG.level_tiles[dim]
    print('number of tiles', DZG_dim)
    if debug:
        print('DZGdim: ', DZG_dim)
    if debug:
        print('pixels:',pixels)
        print(pixels[0]/tile_size)

        print(filepath.rstrip(data_format))
    save_directory = save_dir + str(tile_size) + '/'
    directory = '{0}{1}_{2}'.format(save_directory, 
                                    filename,
                                    str(tile_size))

#     if save_tiles: 
#         if not os.path.exists(directory):
#             os.makedirs(directory)

    if debug:
        print(directory)

    # create matrix for writing average values
    avg_matrix = np.zeros((DZG_dim + (3,)))

    for i in range(0,DZG_dim[0]):#pixels[0],tile_size):
        for j in range(0,DZG_dim[1]):#pixels[1], tile_size):
            if debug:
                print(i, j)

            #filename = os.path.basename(file).rstrip('.svs')
            tilename = '{0}/{1}_{2}_{3}.jpg'.format(directory,filename,i,j)
            
            tile = DZG.get_tile(dim,(i,j))
            tile_size_x, tile_size_y = DZG.get_tile_dimensions(dim,(i,j))
            tile_avg = tile.resize((1, 1))
            color = tile_avg.getpixel((0, 0))

#             avg_matrix[i,j,0]=color[0]
#             avg_matrix[i,j,1]=color[1]
#             avg_matrix[i,j,2]=color[2]   
            avg_matrix[i,j]=color
            
#             #print('RGB', R, G, B)
            if (tile_size_x != tile_size) or (tile_size_y != tile_size):
                if debug:
                    print('new_tile_size')

                new_tile = Image.new("RGB", (tile_size, tile_size), color='white')
                new_tile.paste(tile, (0,0))
                tile = new_tile
                if debug: 
                    print(tile.size)
    
    print(avg_matrix)   
    if save_tiles:
        x,y = range(0,DZG_dim[0]), range(0,DZG_dim[1])
        sums = avg_matrix.sum(axis=2)
        maximum = np.max(sums)
        print('sums of colors', sums)
        print('maximum: ', maximum, 'minimum: ', np.min(sums))
        count=0
        count_all=0
        for i in x:
            for j in y:
                count_all+=1
                if not os.path.exists(directory):
                    os.makedirs(directory)
                filename_ = os.path.basename(filename).rstrip(data_format)
                tilename = '{0}/{1}_{2}_{3}.jpg'.format(directory,filename_,i,j)
                if save_empty_tiles:
                    tile.save(tilename, "JPEG")
                else:
                    if sums[i,j] < (threshold*maximum): # changed from 0.95 to 0.97
                        tile.save(tilename, "JPEG")
                        count+=1
        print('number of all tiles: ', count_all, '\n number of saved tiles: ', count)
    slide.close()
    if debug:
        print('avg matrix: ', avg_matrix)
        
        
    if save_tile_map:
        m_new = SaveTileMap(filename, save_dir, avg_matrix, threshold)
        if print_tile_map:
            PrintTileMap(m_new)
    return avg_matrix


def SaveTileMap(filename, path, tilemap, threshold=0.95):

    textfile_path = os.path.join(path, "txt/")

    filename = filename + "_224.txt"
    # filename = '{0}txt/{1}'.format(path,filename)
    if not os.path.exists(textfile_path):
        os.makedirs(textfile_path)
    filename = os.path.join(textfile_path, filename)
    print(filename)
    
    
    x,y = tilemap.shape[0:2]
    m_new = np.zeros((x,y))
    sums = tilemap.sum(axis=2)
    maximum = np.max(sums)

    for i in range(0,x):
        for j in range(0,y):
            if sums[i,j] < (threshold*maximum):
                m_new[i,j] = 1
    np.savetxt(filename, m_new, fmt='%d')
    print('saved tilemap ', filename)
    return m_new

def PrintTileMap(tile_map):
    print(tile_map)
    x_lim, y_lim = tile_map.shape
    # tile_map_size = x_lim, y_lim
    print('tile_map_size', x_lim, y_lim)
    fig1 = plt.figure(figsize=(20,5))
    a=fig1.add_subplot(1,2,1)
    plt.imshow(tile_map.T, cmap='gray_r')
    plt.show()

here is my Julia adaptation:

function image_tiler(filename, filepath, save_dir; 
    tile_size = 224, overlap = 0,save_tiles = false, 
    debug = false, save_empty_tiles = false, 
    save_tile_map = true, threshold = 0.95, print_tile_map = false, data_format = ".mrxs")

    """Returns all tiles of a whole slide image"""
    slide = ops.OpenSlide(filepath)       

    DZG = dz.DeepZoomGenerator(slide, tile_size, overlap, limit_bounds = false)
    # choose a level 
    dim = pyconvert(Int, DZG.level_count)
    # number of pixels in the current dimension(level)
    level_tiles = pyconvert(Tuple, DZG.level_tiles)
    level_dims = pyconvert(Tuple, DZG.level_dimensions)

    pixels = level_dims[dim]

    # number of tiles in the current dimension(level)
    DZG_dim = level_tiles[dim]

    if debug
        println("Slide : $slide \n   File: $filename")
        println("The number of tiles : $DZG_dim" )
        println("The number of pixels : $pixels")
        println("pixel to tile ratio : ", pixels[0] / tile_size)
        println("The file path : ")
    end

    save_directory= "$save_dir" *  string(tile_size) * "/"
    directory = "$save_directory" * "$filename" * "_" * string(tile_size)

    if debug
        println("The save directory : $directory")
    end


    x_dim = DZG_dim[1]
    y_dim = DZG_dim[2]
    # create matrix for writing average values
    avg_matrix = zeros(tuple(DZG_dim..., (3,)...))

    for i in 1:x_dim
        for j in 1:y_dim
            tilename = "$directory" * "/" * "$filename" * "_" * string(i) * "_" * string(j) * ".jpg"
            tile = DZG.get_tile(dim-1, (i-1,j-1))
            tile_size_x, tile_size_y = DZG.get_tile_dimensions(dim-1,(i-1,j-1))
            tile_size_x , tile_size_y =  pyconvert(Int, tile_size_x), pyconvert(Int, tile_size_y)

            tile_avg = tile.resize((1,1))
            color = tile_avg.getpixel((0,0))
            avg_matrix[i,j,:] = pyconvert(Array, color)

            if (tile_size_x != tile_size) || (tile_size_y != tile_size)

            new_tile = Image.new("RGB", (100,100), color="white")
            new_tile.paste(tile, (0,0))
            tile = new_tile
                if debug
                    println(tile.size)
                end
            end
        end
    end

    if save_tiles
        sums = sum(avg_matrix, dims = 3)
        max_sum = maximum(sums)
        println("Sums of colors: $sums")
        println("Max sum: $max_sum    Min sum: $(minimum(sums))")
        count = 0
        count_all = 0
        for i in 1:x_dim
            for j in 1:y_dim
                if !(isdir(directory))
                    mkpath(directory)
                end
                filename = basename(filename) |> splitext
                filename = filename[1]
                tilename = "$directory" * "/" * "$filename" * "_" * string(i) * "_" * string(j) * ".jpg"
                if save_empty_tiles
                    tile.save(tilename, "JPEG")
                else
                    if (sums[i,j] > threshold * max_sum)
                        tile.save(tilename, "JPEG")
                        count += 1
                    end
                end
                    
         
                count_all += 1            
            end
        end

    end

    slide.close()

    avg_matrix = avg_matrix/255
    if save_tile_map
        m_new = save_tile_map(filename, directory, avg_matrix, threshold)
        if print_tile_map
            display_tilemap(m_new)
        end
    end

    return avg_matrix
end

function save_tile_map(filename, path, tilemap, threshold=0.95)

    textfile_path = joinpath(path, "txt/")
    filename = filename * "_224.txt"

    #create text file directory
    if !(isdir(textfile_path))
        mkdir(textfile_path)
    end

    filename = joinpath(textfile_path, filename)
    println(filename)

    sums = reshape(sum(tilemap,dims=3), (size(tilemap)[1], size(tilemap)[2]))
    max_sum = maximum(sums)

    m_new = sums .< threshold * max_sum

    writedlm(filename, m_new)
    println("Saved tile map: $filename    with threshold: $threshold")

    return m_new
end

function display_tilemap(tile_map; color=true)
    imshow(M) = get.(Ref(ColorSchemes.rainbow), M ./ maximum(M))
	imshow(x::AbstractVector) = imshow(x')
    # println(tile_map)
    x_lim, y_lim = size(tile_map)
    println("tile map size :  $x_lim   $y_lim")
    if color 
        imshow(tile_map)
    else
        Gray.(tile_map)
    end
  
end
2

If I’m not mistaken, it looks like you’re primarily calling some existing python code from your Julia implementation? In that case you’re probably mostly seeing the overhead of communication between the two languages.

2 Likes
3

Thank you! I was thinking it may speed up at least the for loop

4

Ah yes, good idea but for that to work you’d probably have to have no python code inside the loop.

1 Like
5

In particular, you have to be careful about your datastructures to ensure that no copies are made going to/from Python. pyconvert(Array, color) makes a copy, for example.

But yes, you need to remove Python from your innermost loop to have any hope of a speedup.

3 Likes
6

Thank you! Unfortunately Julia does not have openslide, so I have to use python code inside the loop.

7

Thank you for your kind suggestion! If I do not use pyconvert to convert python to Julia, what else could I do?

8

Those loops don’t look particularly hot, so I doubt there’s all that much time spent looping in the Python code. I would be surprised if you can get any speedup at all without rewriting significant parts of the Python functions you’re calling as Julia code, but you should really profile the code to see where the time is actually spent. And if most of the time is spent in C code or similar, you likely need to come up with algorithmic improvements or good use of multithreading and/or SIMD for a chance to improve.

Depends on whether you’re comfortable with ccall-ing the OpenSlide library yourself.

julia> using OpenSlide_jll

julia> slide = ccall((:openslide_open, libopenslide), Ptr{Cvoid}, (Cstring,), "x.mrxs")
Ptr{Nothing} @0x00000000019cc000

julia> width = Ref{Int64}(0)
Base.RefValue{Int64}(0)

julia> height = Ref{Int64}(0)
Base.RefValue{Int64}(0)

julia> ccall((:openslide_get_level_dimensions, libopenslide), Cvoid,
             (Ptr{Cvoid}, Int32, Ref{Int64}, Ref{Int64}),
             slide, 2, width, height)

julia> (width[], height[])
(22848, 50240)

But as already said, if most of the time is spent inside the OpenSlide library, it won’t matter whether you call it from Python or from Julia.

3 Likes
9

Thank you for your informative reply! I found almost all of the time spend in the following part. I am not sure what the ccall is doing. Appreciate your help though! :+1:

10

The inside of that loop has many different concerns - something that is repeated throughout your code

if you wrote pseudocode of the operations

for every i
   for every j
       create a tilename
       get the tile
       determine the size
       record the average color
       conditionally replace the current tile
   end j
end i

if you then break out each of those concerns into functions, you can reason about each one at a time and, more importantly, see which one takes the most time.

as a side note,

tilename = joinpath(directory, "$(filename)_$(i)_$(j).jpg")

would be better, but that contributes very little overhead to a 1hr runtime! :slight_smile: - it is also not used

As far as I can see, the current tile is discarded on every iteration of the loop.

This approach also reveals that there is no data dependence between iterations, so one could possibly parallelize that - although that assumes the python stuff is threadsafe ( [ X ] to doubt ).

At the very least, one could spawn a new task to calculate the average color, and continue on reading the next tile.

Something like this (untested)

function avgcolor(tile)
    tile_avg = tile.resize((1,1))
    color = tile_avg.getpixel((0,0))
    pyconvert(Array, color)
end

@sync for i in 1:x_dim
    for j in 1:y_dim
        ...
        @async avg_matrix[i,j,:] = avgcolor(tile)
        ...
    end
end

Anyway, if you separate stuff out, you can get a better idea of what needs improving and have a vision of how to and even if one can move forward

2 Likes
11

Thank you so much for your very informative reply! :smiley:

I will try that out!

basically, What I expected it would do is to patch the tile when it is not the same as others, not discarding it.

1 Like