Title: Tutorial 1: High-level programming in Quasar

Convolution: from Matlab to Quasar

In this example we will port a simple Matlab function to a Quasar function. Lets start from following Matlab code:

% load image and cast to double precision
im = double(imread('lena_big.tif'));

% initialize variables
y = zeros(size(im));
s = 0;
channel = 2;

% define filter mask
N = 5;
mask = ones(2*N+1,2*N+1)/(2*N+1)^2;
dim = size(im);

% filtering through straightforward convolution implementation and
% summation over specified channel
tic();
for m=1:dim(1)
    for n=1:dim(2)
        % summation
        s = s + im(m,n,channel);
        
        % filter only over inner pixels to avoid boundary issues
        if (n>N && m>N && n<size(im,1)-N+1 && m<size(im,2)-N+1)
            a = zeros(1,1,3);
            for k=-N:N
                for l=-N:N
                    a = a + mask(1+N+k,1+N+l) * im(m-k,n-l,1:3);
                end
            end
            y(m,n,1:3) = a;
        end
    end
end
toc();

% show results
disp(num2str(s));
disp(num2str(sum(sum(im(:,:,channel)))))
imshow(uint8(y));

The above code applies a linear filter to an image through a straightforward convolution and (meanwhile) computes the sum of intensities in a specified channel. In the filtered image, a pixel value corresponds to a weighted average of its surrounding pixels in the input image. In order to port this small script to Quasar we can practically copy the code, with some minor adjustments:

1. Matrix/vector indices start with 0 rather than with 1.
2. Idem for built-in functions that require dimension indices (e.g., size(im,1) becomes size(im,0)).
3. A matrix element is indicated using square brackets instead of round, i.e. M[i,j] instead of M(i,j)
4. Strings are indicated by double quotes vs. single quotes in Matlab
5. For a more elaborate comparison between Matlab and Quasar we refer to: Comparison

There are some motivations for these syntax choices in Quasar: first, indices start with 0 because of performance reasons (especially when generating code for the GPU; in some cases it is difficult to optimize the base-1 index offseting away). Second, matrices use brackets in Quasar to make the difference between function calls and matrix/vector indexing more clear.

These adjustments result in the following Quasar code:

function [] = main()
    % load image and cast to double precision
    im = imread("lena_big.tif")

    % initialize variables
    y = zeros(size(im))
    s = 0
    channel = 1

    % define filter mask
    N = 5
    mask = ones(2*N+1,2*N+1)/(2*N+1)^2
    dim = size(im)

    % filtering through straightforward convolution implementation and
    % summation over specified channel
    tic();
    for m=0..dim[0]-1
        for n=0..dim[1]-1
            % summation
            s += im[m,n,channel]
            
            % filter only over inner pixels to avoid boundary issues
            if n>=N && m>=N && n<size(im,0)-N && m<size(im,1)-N
                a = [0.0,0.0,0.0]
                for k=-N..N
                    for l=-N..N
                        a += mask[N+k,N+l] * im[m-k,n-l,0..2]
                    end
                end
                y[m,n,0..2] = a
            endif
        end
    end
    toc()

    % show results
    print(s)
    print(sum(im[:,:,channel]))
    imshow(y)
end

Implementation using a kernel function

The code can also be rewritten such that the for loops than iterate over all pixels are replaced by a kernel function that runs on each pixel. On its own, this is not a necessity (because the automatic loop parallelizer performs exactly this task!), but it demonstrates how kernel functions that are typically programmed in CUDA/OpenCL can directly be written from within the host language. This results in the following Quasar code:

function [] = main()
    image = imread("lena_big.tif")/255
    
    r = 5
    filt = ones(2*r+1,2*r+1)/(2*r+1)^2
    
    [M,N] = size(image,0..1)
    
    channel = 2
    
    out = uninit(size(image))
    
    function [s:scalar] = __kernel__ filter_kernel(image:cube, filt:mat, pos:ivec2)
        s += image[pos[0],pos[1],channel]
        
        cumulative = [0.0,0.0,0.0]
        for k = -r..r
            for l = -r..r
                cumulative = cumulative + filt[k-r,l-r] * image[pos[0]+k,pos[1]+l,0..2]
            end
        end
        
        out[pos[0],pos[1],0..2] = cumulative
    end
    tic()
    s = parallel_do([M,N],image,filt,filter_kernel)
    toc()
    
    print(s)
    print(sum(image[:,:,channel]))
    imshow(out)
end

Boundary handling

The Matlab script avoids going out of boundaries by only calculating the filter for the interior of the image, and by ignoring the border of the image. However, for many filtering operations it is assumed that outside the image all values are zero, or a mirroring of the values along the boundary, or that an image is periodic (i.e., im[-1,1]=im[M-1,1], with M the width of the image). Normally these types of boundary handling would require additional coding. Within Quasar these extension can be easily achieved using the right modifiers. Adjust the auto_parallel script such that the filter is also calculated near the boundary. Change the boundary handling, e.g. switch from zeros to mirroring

function [s:scalar] = __kernel__ filter_kernel(image:cube'mirror, filter:mat, pos:ivec2)
    s += image[pos[0],pos[1],channel]
    
    cumulative = [0.0,0.0,0.0]
    for k = -r..r
        for l = -r..r
            cumulative = cumulative + filter[k+r,l+r] * image[pos[0]+k,pos[1]+l,0..2]
        end
    end
    
    out[pos[0],pos[1],0..2] = cumulative
end