Adding identification spots to a sigma clipped image

I have come across the following from IA in which he demonstrates how to perform sigma clipping for a set of images. I was wondering how to show on the composite image by plotting red circles all those pixels that have been modified. I have tried:

plot(pixelsToReplace,'r.', 'MarkerSize', 6, 'Marker','.');

------------------------------IA's post--------------------------------------------------

The first pass is a for loop to just read in all the images and construct a "mean image" and a "standard deviation image" from the entire set of images. The second pass would be to read in the images again and do the comparison and clipping (if indicated for the pixels in that image). It's just a standard for loop so not sure what the difficulty is. The pseudo code woudl go something like this more or less.... f

or k = 1:numberOfImages
  read in image
  sumImage = sumImage + imageArray
  sum2Image = sum2Image + imageArray .^2;
end
meanImage = sumImage/numberOfImages;
varianceImage = (sum2Image/numberOfImages) - meanImage.^2;
stdDevImage = sqrt(varianceImage)

% Now replace the outliers.
for k = 1:numberOfImages
  % read in image
  outputImage = size(imageArray);
  pixelsToReplace = abs(imageArray - meanImage) > (5 * stdDevImage)
  outputImage(pixelsToReplace) = meanImage(pixelsToReplace)

    %I then add each modified image to a stack called imStack2
end

    %Now create a median of the resultant stack
    imMedSC=median(imStack2,3);
    figure
    imshow(imMedSC,[])

[Adithya Addanki]

Hi Jason,

You may be able to calculate the Image's pixel position using the function "imref2d"/"imref3d". Please refer to the link below:

http://www.mathworks.com/help/images/image-coordinate-systems.html

After fetching the positions of the pixels you may have modified, you may plot them using the "plot" function.

For instance, assuming X and Y are the arrays consisting x and y coordinates of the pixels modified.

% display image using imshow
hold on;
plot(X, Y, 'ro')

These X and Y may need to be calculated using the XWorldLimits and YWorldLimits values returned from imref2d/3d.

Thank you, Adithya