function im = cropxcorr(im_crop, im_orig)
    % CROPXCORR Scaled normalized two-dimensional cross-correlation
    
    % C = CROPXCORR(TEMPLATE,A) computes the normalized cross-correlation of
    % matrices TEMPLATE and A. The matrix A can be smalled than TEMPLATE
    % Returned value C is the calculated TEMPLATE
    
    im_orig_grey = rgb2gray(im_orig);
    [im_orig_grey_h, im_orig_grey_w] = size(im_orig_grey);


    im_crop_grey = rgb2gray(im_crop);
    [im_crop_grey_h, im_crop_grey_w] = size(im_crop_grey);


    %order to speed up processing, scale the maximum image size to 320x240
    %find the constraining dimension
    im_orig_grey_scaler = max(240/im_orig_grey_h,320/im_orig_grey_w);
    im_orig_grey_scaled = imresize(im_orig_grey, im_orig_grey_scaler, 'nearest');
    [im_orig_grey_scaled_h, im_orig_grey_scaled_w, ~] = size(im_orig_grey_scaled);

    %work out which is the constraining dimension
    if im_crop_grey_h/im_orig_grey_scaled_h > im_crop_grey_w/im_orig_grey_scaled_w
        im_orig_grey_dim = im_orig_grey_scaled_h;
        im_crop_grey_dim = im_crop_grey_h;
        accuracy_offset = 6;
        start_offset = 17; % 3*accuracy_offset - 1
    else
        im_orig_grey_dim = im_orig_grey_scaled_w;
        im_crop_grey_dim = im_crop_grey_w;
        accuracy_offset = 8;
        start_offset = 23; % 3*accuracy_offset - 1
    end

    im_dim_max_c = 0;
    im_dim = 0;

    % scale cropped image in increments from 'accuracy_offset' to original image width
    % finding the maximum peak for each new scaled image
    for  i = start_offset:accuracy_offset:im_orig_grey_dim

        c = normxcorr2(imresize(im_crop_grey, i/im_crop_grey_dim, 'nearest'), im_orig_grey_scaled);
        % offset found by correlation
        max_c_temp = max(abs(c(:)));

        if max_c_temp > im_dim_max_c
            im_dim_max_c = max_c_temp;
            im_dim = i;
        end 
    end

    % get the scaler of the maximum
    im_crop_grey_scaler = im_dim/im_crop_grey_dim;

    % compensate by the resize of original image
    im_crop_grey_scaler = im_crop_grey_scaler/im_orig_grey_scaler;


    % we know the rough scale, lets compensate the original image scale
    % to possibly obtain better result

    im_crop_grey_scaled = imresize(im_crop_grey, im_crop_grey_scaler);
    [im_crop_grey_scaled_h, im_crop_grey_scaled_w] = size(im_crop_grey_scaled);

    c = normxcorr2(im_crop_grey_scaled, im_orig_grey);

    % offset found by correlation
    [~, imax] = max(abs(c(:)));
    [ypeak, xpeak] = ind2sub(size(c),imax);
    xoffset = xpeak - im_crop_grey_scaled_w;
    yoffset = ypeak - im_crop_grey_scaled_h;

    
    xoffset  = max(xoffset, 0);
    xbegin = xoffset + 1;
    xend   = xoffset + im_crop_grey_scaled_w;
    yoffset = max(yoffset, 0);
    ybegin = yoffset + 1;
    yend   = yoffset + im_crop_grey_scaled_h;

    if xend > im_orig_grey_w
        xbegin = xbegin - (xend - im_orig_grey_w);
        xend = im_orig_grey_w;
    end
    
    if yend > im_orig_grey_h
        ybegin = ybegin - (yend - im_orig_grey_h);
        yend = im_orig_grey_h;
    end
    
    % extract region
    im = im_orig(ybegin:yend,xbegin:xend,:);
    im = imresize(im, [im_crop_grey_h, im_crop_grey_w]);
end