Equation for 3D Noise Power Spectrum (NPS)

Question

I have this equation for calculating 3D NPS. I need to obtain to obtain 2D planar NPS by integrating the 3D NPS along the the frequency axis orthogonal to the frequency plane of interest. I developed this code using matlab as below. I have also attached the equation I should use to calculate the 3D NPS as an image file. Can someone please help in identifying what is wrong with this code? My error message is : integral_value = integral3(integrand, x_lim(1), x_lim(2), y_lim(1), y_lim(2), z_lim(1), z_lim(2), 'Method', 'iterated');

%voxel pitch
    DeltaX = 0.1; 
    DeltaY = 0.1; 
    DeltaZ = 0.1; 

    %Nx = 400;
    Nx = 400;
    Ny = 400; 
    Nz = 2;    %this represents the number of voxels in each dim

  %  NumImages = 50;
 %   voiSize = 40;
   

%Define a box volume around the centroid and finding the indices of the box
centroid = [637, 641];
%halfSize = floor(voiSize / 2);
voxeldim = [400,400];
%halfSize = floor(VoiSize/ 2);
 halfSize = floor(voxeldim/2);
voiBoundingBox = [
    centroid(1) - halfSize(1) + 1, centroid(1) + halfSize(1);
    centroid(2) - halfSize(2) + 1, centroid(2) + halfSize(2)
];

% Initialize NPS
 
NPSxyz = zeros(Nx,Ny,Nz);

meanVOIValue = 57.8085           % for size of 30 slices
%d_value = zeros(400,400);
d_value = zeros(400,400);

fx = (-Nx/2:Nx/2-1) / (Nx * DeltaX);  %where fx,fy and fz are the spatial frequencies in x,y and z axis  (are these the sampling frequencies along x,y and z axis)?
fy = (-Ny/2:Ny/2-1) / (Ny * DeltaY);
fz = (-Nz/2:Nz/2-1) / (Nz * DeltaZ);



dirName = pwd;              %# folder path
files = dir( fullfile(dirName,'*') );  
files = {files.name}';                      %'# file names
files = files(3:end)

 integral_result = 0;
%integral_result = 0;
%Loop over the stack of images
for i = 1:Nz
     con= 520+i;
    imgint(:,:,:,i)=squeeze(mat2gray(dicomread(fullfile(dirName,files{con}))));
    img=imgint;
    %slimg = img(:,:,3,:);
    %image = imshow(img(:,:,i));
    NPSimg = img(:,:,1,i);
    %convert image to 32 bits
    int32Data = int32(NPSimg);

    %Extract the volume
    voiImage = int32Data(voiBoundingBox(1, 1):voiBoundingBox(1, 2), ...
                      voiBoundingBox(2, 1):voiBoundingBox(2, 2));
   
    %FT from spatial to freq. domain for the extracted volume
     fftData = fftshift(fftn(voiImage));
     d_value = abs(fftData);
       % Mean Vol value across the stack
     
    %meanVoxelValues(i) = mean(abs(fftData(:)));

    % Create frequency grid
    %[Nx, Ny, Nz] = size(NPSimg);
    integral_result = 0;        

%NPSxyz = zeros (Nx, Ny, Nz);
%for z = 1:2 %Nz
    for x = 1:Nx
        for y = 1:Ny
            
            mean_d = meanVOIValue;
            integrand = @(x, y, z)(d_value - mean_d) .* exp(-2i * (fx(x) .* x + fy(y) .* y + fz(i) .* i));
            
            [X, Y] = meshgrid(voiBoundingBox(1, 1):voiBoundingBox(1, 2), ...
                      voiBoundingBox(2, 1):voiBoundingBox(2, 2));
            x_lim = [min(X(:)), max(X(:))];
            y_lim = [min(Y(:)), max(Y(:))];
           % z_lim = [1, size(voiImage, 2)];
             z_lim = [1, Nz];
            % Perform the triple integral
            integral_value = integral3(integrand, x_lim(1), x_lim(2), y_lim(1), y_lim(2), z_lim(1), z_lim(2), 'Method', 'iterated');
             modulated_result = abs(integral_value).^2;
             % Normalize by voxel pitch and volume size
             NPSxyz(:, :,i) = (DeltaX * DeltaY * DeltaZ) / (Nx * Ny * Nz) * modulated_result;

            
        end
    end
   



    
  % NPSxy(:, :,i) = (DeltaX * DeltaY)/(Nx*Ny)* modulated_result;
end

% Integrate along the frequency axis orthogonal to the z-axis
    integrated_result = sum(NPSxyz, 3);

figure 
% Display the logarithm of the NPS for better visualization
new_img = log10(1+integrated_result);
imagesc(new_img);
%imagesc(fx,fy,integrated_result);
axis xy;    % Set the y-axis to be increasing from bottom to top
%imagesc(integrated_result);

title('2D Noise Power Spectrum','FontName', 'Helvetica', 'FontSize', 14);
xlabel('Spatial Frequency (cycles/mm)','FontName', 'Helvetica', 'FontSize', 13);
ylabel('Spatial Frequency (cycles/mm)', 'FontSize', 13);
colorbar;
colormap('parula');