function output = generate_simulated_spots(npixels, digitizeFLAG, ... IntensityNoiseFLAG, positionNoiseFLAG) % % a function to generate an image of spots to centroid and to provide the original % locations of the spots. % % input % ----- % npixels - pixels in output image % digitizeFLAG: 0 - floating point array is output % 1 - intensity array is rounded down to nearest % integer (12-bit) % IntensityNoiseFLAG: 0 - no photo-electron shot noise % 1 - photo-electron shot noise is added % positionNoiseFLAG: 0 - position noise off % : 1 - position noise ON - deviation = positionDeviationP % output % ------ % output.image - the 2D intensity array [npixels, npixels] % output.positions - the nominal positions of the spots: [N, 2] array % % EXAMPLE % -------- % % output = generate_simulated_spots(1024, 1, 1, 0, 0) % - generates an image with photo-electron intensity noise, but with % the theoretical spot positions undeviated from a regular array. Fast % Gaussian intensity calculated. % % Gaussian Integration % --------------------- % % - the intensity, Ip, in a pixel is given by the following % Ip = integrate( exp{ -2 [(x-xc)^2 + (y-yc)^2]/w^2 }, {x, xa, xb}, {y, ya, yb} ) % = 1/8 * pi * w^2 % * [ erf( sqrt(2)*(x0-xc)/w ) - erf( sqrt(2)*(x1-xc)/w ) ] % * [ erf( sqrt(2)*(y0-yc)/w ) - erf( sqrt(2)*(y1-yc)/w ) ] % % % Coordinate representation of pixels % ----------------------------------- % % | | | % | | | % 2 +----------+----------+--- % | | | % | | | % | | <========== this pixel has indicies [2, 2] but is % | | | said to have coordinates [1.5, 1.5], % | | | corresponding to the center. % 1 +----------+----------+--- % | | | % | | | % | | | % | | | % | | | % 0 +----------+----------+--- % 0 1 2 % % % Written by Aidan Brooks. 19th Mar 2010 % define the physical parameters of the HWS pixelSize = 12.0E-6; % pixel size in m maxElectronsPerPixel = 1.6E5; % around 160,000 electrons before saturation maxElectronsInExposure = 1.5E5; holePitchM = 430E-6; % separation of the holes in the Hartmann plate holePitchP = holePitchM/pixelSize; % separation of the holes in the Hartmann plate in Pixels spotRadiusM = 75.0E-6; % approximate spot radius on CCD = 1/2 hole size spotRadiusP = spotRadiusM/pixelSize; CCDbackground = 50.0; % background level in 12-bit representation positionDeviationP = 0.1; % the standard deviation of random spot position deviation (in pixels) intensityUnitsPerElectron = 4096/maxElectronsPerPixel; % buffer (in pixels) around the edges buffer = ceil(holePitchP); % generate the array of positions and the image if the array is big enough if (npixels - 2*buffer - holePitchP > 0) % generate the blank image image_arr = zeros(npixels, npixels); % generate the coordinate arrays - the pixel value is the center of the % pixel, therefore the first pixel is 0.5, the next 1.5, etc. x = 1:npixels; x2D = (x'*ones(1,npixels))' - 0.5; y2D = x2D'; % get the number of spots across the top of a recti-linear array nspot1D = floor((npixels - 2*buffer)/holePitchP) + 1; % create an empty array of the nominal spot positions nominal_positions = []; % loop through the spots for ii = 1:nspot1D % ii-th column for jj = 1:nspot1D % jj-th row spot_position = ([(ii - 1), (jj-1)]*holePitchP) + [buffer buffer]; % if its flag is present then add position noise if (positionNoiseFLAG ~= 0) spot_position = spot_position + ... random('Normal', 0, positionDeviationP, 1, 2); end % add the spot_position to the nominal position list nominal_positions = [nominal_positions; spot_position]; %----------------------------------------------- % IMAGE % in a square around the hole, add the intensity xc = spot_position(1); yc = spot_position(2); imin = floor(xc - holePitchP/2); imax = ceil(xc + holePitchP/2); jmin = floor(yc - holePitchP/2); jmax = ceil(yc + holePitchP/2); % get the limits of integration % NOTE: [row number (= Y), column number (= X)] x0 = x2D(jmin:jmax, imin:imax) - 0.5; x1 = x2D(jmin:jmax, imin:imax) + 0.5; y0 = y2D(jmin:jmax, imin:imax) - 0.5; y1 = y2D(jmin:jmax, imin:imax) + 0.5; % get the integrated gaussian using erf intGauss = (1/8)*pi*spotRadiusP^2 * ... ( erf( sqrt(2) * (x0 - xc)/spotRadiusP ) - ... erf( sqrt(2) * (x1 - xc)/spotRadiusP )).* ... ( erf( sqrt(2) * (y0 - yc)/spotRadiusP ) - ... erf( sqrt(2) * (y1 - yc)/spotRadiusP )); % add the integrated sub-array to the total image image_arr(jmin:jmax, imin:imax) = image_arr(jmin:jmax, imin:imax) + intGauss; %----------------------------------------------- end end % rescale the image to be in photoelectrons image_arr = maxElectronsInExposure*image_arr/max(max(image_arr)); % in photo-electrons % add the intensity noise - a normal distribution where the deviation in photo-electrons is % proportional to the sqrt of the number of the photo-electrons if (IntensityNoiseFLAG ~= 0) image_arr = image_arr + ... random('Normal', 0, sqrt(image_arr), npixels, npixels); end % rescale the image to be in 12-bit digital units image_arr = image_arr*intensityUnitsPerElectron; % add the background level image_arr = image_arr + CCDbackground; % digitize the image if this is requested. if (digitizeFLAG ~= 0) image_arr = floor(image_arr); end output.positions = nominal_positions; output.image = image_arr; else disp(['The number of pixels must be greater than ', num2str(ceil(2*buffer + holePitchP))]) output = 0; end