How to calculate noise power spectrum of an noise image

I have an uniform image 512 512 pixels , I want to calculate a ROI(128 128 pixels) in this image by using this formula:

u and v are spatial frequency (mm−1) in the x and y directions, respectively, dx and dy are pixel size (mm), Nx and Ny are the number of pixels in the x and y direction of the ROI, F[] denotes the 2D Fourier transform, I(x,y) is the pixel value (HU) of a ROI at position (x,y), and P(x,y) is a 2nd order polynomial fit of I(x,y).

My questions is: 1/ how can I make a ROI in the center of image with size 128 128 2/ how can I take values of pixels in ROI for using like I(x,y) 3/ how can I evaluate 2nd order polynomial fir of I(x,y), and 3/ how can I calculate NPS with all of these

Thanks you for your helping.

Answer 1

[Image Analyst]

To get the middle pixel you can do

[rows, columns, numberOfColorChannels] = size(yourImage);
middleRow = floor(rows/2);
middleColumn = floor(columns/2);

To crop out +/- N pixels around that do

leftColumn = middleColumn - N;
topRow = middleRow - N;
croppedImage = imcrop(yourImage, [leftColumn, topRow, 2*N, 2*N]);

With a uniform image your FFT will be a very narrow sinc function.

[Phuong Phan Hoai]

Thanks . But If I want to crop a region with size 64x64 or 128x128 pixels in every where on an image, how can I do that?

[Image Analyst]

I don't know what "in every where on an image" means.

imcrop takes the bounding box in the form [leftColumn, topRow, width, height].

You can make any of those whatever you want to locate the box in the desired location.

Answer 2

[Sang Hyeok Park]

The following is the Matlab code for noise power spectrum.

%Clearing the memory and screen values clc;clear;

% Reading in dicom flat field image

info = dicominfo('I0001_1'); A = dicomread(info); %figure(1);imagesc(A); 91

%Splitting the image into 128 x 128 regions and taking the Fourier Transform of each section

F=zeros(128,128);

for i=500:128:1396

for j=500:128:1396

T = A(i:i+127,j:j+127);

P = (log(abs(fft2(T))).^2);

F = F+P; end end

code from http://purl.fcla.edu/fcla/etd/UFE0001030

I think this code may be helpful to your question.