Here's an alternative without convolution. One intricacy that you must take into account is that if you want to place each image at the centre of each dot, you must determine where the top left corner is and index into your output image so that you draw the desired object from the top left corner to the bottom right corner. You can do this by taking each black dot location and subtracting by half the width horizontally and half the height vertically.

Now onto your actual problem. It's much more efficient if you loop through the set of points that are black, not the entire image. You can do this by using the find command to determine the row and column locations that are 0. Once you do this, loop through each pair of row and column coordinates, do the subtraction of the coordinates and then place it on the output image.

I will impose an additional requirement where the objects may overlap. To accommodate for this, I will accumulate pixels, then find the average of the non-zero locations.

Your code modified to accommodate for this is as follows. Take note that because you are using JPEG compression, you will have compression artifacts so regions that are 0 may not necessarily be 0. I will threshold with an intensity of 128 to ensure that zero regions are actually zero. You will also have the situation where objects may go outside the boundaries of the image. Therefore to accommodate for this, pad the image sufficiently with twice of half the width horizontally and twice of half the height vertically then crop it after you're done placing the objects.

mainImage=imread('https://i.stack.imgur.com/gbhWJ.png');
secondaryImage=imread('https://i.stack.imgur.com/P0meM.png');
secondaryImageResized = imresize(secondaryImage, [250 300]);

% Find half height and width
rows = size(secondaryImageResized, 1);
cols = size(secondaryImageResized, 2);
halfHeight = floor(rows / 2);
halfWidth = floor(cols / 2);

% Create a padded image that contains our main image.  Pad with white
% pixels.
rowsMain = size(mainImage, 1);
colsMain = size(mainImage, 2);
outputImage = 255*ones([2*halfHeight + rowsMain, 2*halfWidth + colsMain, size(mainImage, 3)], class(mainImage));
outputImage(halfHeight + 1 : halfHeight + rowsMain, ...
        halfWidth + 1 : halfWidth + colsMain, :) = mainImage;

% Find a mask of the black pixels
mask = outputImage(:,:,1) < 128;

% Obtain black pixel locations
[row, col] = find(mask);

% Reset the output image so that they're all zeros now.  We use this
% to output our final image.  Also cast to ensure accumulation is proper.
outputImage(:) = 0;
outputImage = double(outputImage);

% Keeps track of how many times each pixel was hit by the object
% This is so that we can find the average at each location.
counts = zeros([size(mask), size(mainImage, 3)]);

% For each row and column location in the image    
for i = 1 : numel(row)
    % Get the row and column locations
    r = row(i); c = col(i);

    % Offset to get the top left corner
    r = r - halfHeight;
    c = c - halfWidth;

    % Place onto final image
    outputImage(r:r+rows-1, c:c+cols-1, :) = outputImage(r:r+rows-1, c:c+cols-1, :) + double(secondaryImageResized);

    % Accumulate the counts
    counts(r:r+rows-1,c:c+cols-1,:) = counts(r:r+rows-1,c:c+cols-1,:) + 1;
end

% Find average - Any values that were not hit, change to white
outputImage = outputImage ./ counts;
outputImage(counts == 0) = 255;
outputImage = uint8(outputImage);

% Now crop and show
outputImage = outputImage(halfHeight + 1 : halfHeight + rowsMain, ...
        halfWidth + 1 : halfWidth + colsMain, :);
close all; imshow(outputImage);

% Write the final output
imwrite(outputImage, 'finalimage.jpg', 'Quality', 100);

We get:

Edit

I wasn't told that your images had transparency. Therefore what you need to do is use imread but ensure that you read in the alpha channel. We then check to see if one exists and if one does, we will ensure that the background of any values with no transparency are set to white. You can do that with the following code. Ensure this gets placed at the very top of your code, replacing the images being loaded in:

mainImage=imread('https://i.stack.imgur.com/gbhWJ.png');

% Change - to accommodate for transparency
[secondaryImage, ~, alpha] = imread('https://i.imgur.com/qYJSzEZ.png');
if ~isempty(alpha)
    m = alpha == 0;
    for i = 1 : size(secondaryImage,3)
        m2 = secondaryImage(:,:,i);
        m2(m) = 255;
        secondaryImage(:,:,i) = m2;
    end
end

secondaryImageResized = imresize(secondaryImage, [250 300]);

% Rest of your code follows...
% ...

The code above has been modified to read in the basketball image. The rest of the code remains the same and we thus get:

You can use convolution to achieve the desired effect. This will place a copy of im everywhere there is a black dot in imz.

% load secondary image
im = double(imread('secondary.jpg'))/255.0;

% create some artificial image with black indicators
imz = ones(500,500,3);
imz(50,50,:) = 0;
imz(400,200,:) = 0;
imz(200,400,:) = 0;

% create output image
imout = zeros(size(imz));
imout(:,:,1) = conv2(1-imz(:,:,1),1-im(:,:,1),'same');
imout(:,:,2) = conv2(1-imz(:,:,2),1-im(:,:,2),'same');
imout(:,:,3) = conv2(1-imz(:,:,3),1-im(:,:,3),'same');
imout = 1-imout;

% output
imshow(imout);

Also, you probably want to avoid saving main.jpg as a .jpg since it results in lossy compression and will likely cause issues with any method that relies on exact pixel values. I would recommend using .png which is lossless and will also likely compress better than .jpg for synthetic images where the same colors repeat many times.