画像にノイズを混ぜる、ぼかす

    gaussSigma = 7;

    filtersize = 2*ceil(2*gaussSigma)+1 % https://jp.mathworks.com/help/images/ref/imgaussfilt.html

    GF = fspecial('gaussian', filtersize, gaussSigma); % Gaussian Filter

    TL = imfilter(TF, GF);

% Original paper

% Seeing Jesus in toast: Neural and behavioral correlates of face pareidolia

% https://www.sciencedirect.com/science/article/pii/S0010945214000288

%

% Please be careful that this code was not written by the authors.

% Daiichiro Kuroki who worked as a technical staff in Kyushu University

% made it by reading the paper and asking the authors.

% The authors and I do not grantee that the program works properly.

%

% Before running the program, face and letter images without noise must be

% saved in two separate folders.

% Be careful that the face-noise and letter-noise images could not be made

% at the same time. You have to set the inputFolder respectively.

%

% You have to change the extension of image files according to your

% stimuli. i.e., imgFileList.

% Moreover, please also check 'stimType' and 'difficulty'.

%

% This program could run in Octave under Ubuntu.

% If you use Matlab, both 'the Statistics and Machine Learning Toolbox' and

% 'the Image Processing Toolbox' are needed.

%

% It takes about 1 minute to make a noise image!

function makeImgWithNoise4

pkg load statistics image

%pkg load statistics

clear all;

% Be careful. Please check also Sigma, lambda, and gamma

inputFolder = ['images' filesep]; % folder name

imgFileList = dir([inputFolder '*.jpg']); % bmp, jpg, jpeg, or png?

stimType = 1; % 1=face, 2=charactor, 3=checker board

difficulty = 1; % 1=easy to detect, 2=hard to detect 

outputFolder = 'output';

if ~exist(outputFolder, 'dir')

    mkdir(outputFolder);

end

checkerBoardWidth = 40; % mod(noiseWidth, checkerBoardWidth) must be zero.

checkerColor1 = [100 100 100];

checkerColor2 = [200 200 200];

% These width/height should be greater than the width/heigth of images.

noiseWidth = 480;

noiseHeight = 480;

%noiseWidth = 800;

%noiseHeight = 700;

% PTB-3 correctly installed and functional? Abort otherwise.

AssertOpenGL;

%Screen('Preference', 'SkipSyncTests', 1); 

% Select screen with maximum id for output window:

screenid = max(Screen('Screens'));

% Open a fullscreen, onscreen window with gray background. Enable 32bpc

% floating point framebuffer via imaging pipeline on it, if this is possible

% on your hardware while alpha-blending is enabled. Otherwise use a 16bpc

% precision framebuffer together with alpha-blending. We need alpha-blending

% here to implement the nice superposition of overlapping gabors. The demo will

% abort if your graphics hardware is not capable of any of this.

PsychImaging('PrepareConfiguration');

PsychImaging('AddTask', 'General', 'FloatingPoint32BitIfPossible');

%[win winRect] = PsychImaging('OpenWindow', screenid, 128);

% Color value seems to be from 0 to 255 not from 0.0 to 1.0 at least in Linux.

% To present a sin-wave gabor, the background color should be grey.

%[win winRect] = PsychImaging('OpenWindow', screenid, 0);

[win, winRect] = PsychImaging('OpenWindow', screenid, 0, [0 0 1000 800] );

% Retrieve size of window in pixels, need it later to make sure that our

% moving gabors don't move out of the visible screen area:

[w, h] = RectSize(winRect);

centerX = w/2;

centerY = h/2;

% Query frame duration: We use it later on to time 'Flips' properly for an

% animation with constant framerate:

ifi = Screen('GetFlipInterval', win);

nImages = size(imgFileList, 1)

if stimType == 3 % checker board

  if mod(noiseWidth, checkerBoardWidth) ~= 0

    error('Please check the checkerBoardWidth');

  end

  if mod(noiseHeight, checkerBoardWidth) ~= 0

    error('Please check the checkerBoardWidth');

  end

  for j = 1:noiseHeight/checkerBoardWidth

    if mod(j,2) == 0

      colorType = 1;

    else

      colorType = 2;

    end

    for i = 1:noiseWidth/checkerBoardWidth   

      tmpRect = [(i-1)*checkerBoardWidth, (j-1)*checkerBoardWidth, i*checkerBoardWidth, j*checkerBoardWidth];

      if colorType == 1

        tmpColor = checkerColor1;

        colorType = 2;

      else

        tmpColor = checkerColor2;

        colorType = 1;

      end

      Screen('FillRect', win, tmpColor, tmpRect);

    end

  end

  Screen('Flip', win);

  imageArray = Screen('GetImage', win, [0, 0, noiseWidth, noiseHeight]);

  imwrite(imageArray, [outputFolder filesep 'checkr-board.jpg']);

  sca;

  return

end

for k = 1:nImages

  imgFileName = [inputFolder imgFileList(k).name]

  x1 = 1:noiseWidth;

  x2 = 1:noiseHeight;

  [X1, X2] = meshgrid(x1, x2);

  vNoise = zeros(length(X1(:)), 1); % Initialize (vNoise is equal to zeros(length(X2(:)), 1)

  for i = 1:3

      % The 'sd' written in the paper seems to be actually covariance.

      if i == 1

        sd = 64;

        contrast = 1;

      elseif i == 2

        sd = 256;

        contrast = 5;

      else

        sd = 1024;

        contrast = 25;

      end

      Sigma = eye(2).*sd;

      ngabors = round((noiseWidth/(sd*0.15))^2);

      %ngabors = 10;

      for j = 1:ngabors

        mu = [rand(1)*noiseWidth, rand(1)*noiseHeight];

        F = contrast.*mvnpdf([X1(:), X2(:)], mu, Sigma);

        vNoise = vNoise + F;

      end

      % for checking

%      noiseMatrix = reshape(vNoise, length(x2), length(x1))./max(vNoise);

%      gabortex=Screen('MakeTexture', win, noiseMatrix, [], [], 2);

%      Screen('DrawTextures', win, gabortex);  

%      Screen('Flip', win);

%      KbStrokeWait;

  end

  noiseMatrix = reshape(vNoise, length(x2), length(x1))./max(vNoise); % normalize from 0 to 1

  NI = (1 - noiseMatrix).*255; % reverse

  imdata = imread(imgFileName);    

  TF = zeros(noiseHeight, noiseWidth); % True Face

  LX = round((noiseWidth - size(imdata,2))/2)+1;

  LY = round((noiseHeight - size(imdata,1))/2)+1;

  %TF(LY:LY+size(imdata,1)-1, LX:LX+size(imdata,2)-1) = imdata(:,:,1);

  TF(LY:LY+size(imdata,1)-1, LX:LX+size(imdata,2)-1) = imdata;

  %imagetex = Screen('MakeTexture', win, TF);

  %Screen('DrawTexture', win, imagetex);

  %Screen('Flip', win);

  %KbStrokeWait;  

  if stimType == 1 % face 

    if difficulty == 1 % easy to detect

      Sigma = eye(2).*2400;

      lambda = 0.3; 

    else % hard to detect

      Sigma = eye(2).*(10^9);

      lambda = 0.9999997; 

    end

    mu = [noiseWidth/2, noiseHeight/2];

    F = mvnpdf([X1(:), X2(:)], mu, Sigma);

    F = F./max(F);

    F = reshape(F, length(x2), length(x1));

    gabor_brob = F;

    whos m gabor_brob NI TF; 

    B = lambda.*(1-gabor_brob)+(1-lambda).*gabor_brob;

    FN = lambda.*NI.*(1-gabor_brob) + (1-lambda).*TF.*gabor_brob;

    FN = FN./B;

    outputData = FN;

  elseif stimType == 2 % charactor

    if difficulty == 1 % easy to detect

      gamma = 0.6;

    else

      gamma = 0.35; % hard to detect (still easy)

    end

    % Blurring a letter. Although this is not written in the paper, the authors told me directly.

    gaussSigma = 7;

    filtersize = 2*ceil(2*gaussSigma)+1 % https://jp.mathworks.com/help/images/ref/imgaussfilt.html

    GF = fspecial('gaussian', filtersize, gaussSigma); % Gaussian Filter

    TL = imfilter(TF, GF);

    %TL = TF; % without blurring

    LN = NI - (158-TL).*gamma; % Negative values are possible.

    outputData = uint8(LN);

  else

  end

  newTex = Screen('MakeTexture', win, outputData);

  Screen('DrawTexture', win, newTex);

  Screen('Flip', win); % image with gabor

  %KbStrokeWait;

  %imageArray = Screen('GetImage', win, [centerX-noiseWidth/2, centerY-noiseHeight/2, centerX+noiseWidth/2, centerY+noiseHeight/2]);

  [pathstr, imgFilename2, ext] = fileparts(imgFileName);

  imwrite(uint8(outputData), [outputFolder filesep imgFilename2 '-withNoise.jpg']);

  Screen('Close');

  %error('stop');

end

% Close onscreen window, release all ressources:

sca;

% Done.

return;

% Little helper function: Converts degrees to radians, included for

% backward compatibility with old Matlab versions:

% function radians = deg2rad(degrees)

%     radians = degrees * pi / 180;

% return