Блог

Я генерирую 3d фрактальный шум в MATLAB, используя различные методы. Это работает относительно хорошо, но у меня проблема, из-за которой я вижу артефакты вертикальной чередования в моем шуме. Это происходит независимо от того, какой тип данных или разрешение я использую.

Редактировать: я понял это. Решение опубликовано в качестве ответа ниже. Спасибо всем за ваши мысли и руководство!

 expo = 2^6;
dims = [expo,expo,expo];
beta = -4.5;

render = randnd(beta, dims); % Create volumetric fractal
render = render - min(render); % Set floor to zero
render = render ./ max(render); % Set ceiling to one
%render = imbinarize(render); % BW Threshold option
render = render .* 255; % For greyscale

slicer = 1; % Turn on image slicer/saver
i = 0; % Page counter
format = '.png';
imagename = '___testDump/slice';

imshow(render(:,:,1),[0 255]); %Single test image

if slicer == 1
    for c = 1:length(render)
    i = i 1;
    pagenumber = num2str(i);
    filename = [imagename, pagenumber, format];
    imwrite(uint8(render(:,:,i)),filename)
    end
end

function X = randnd(beta,varargin)

seed = 999;
rng(seed); % Set seed

%% X = randnd(beta,varargin)

% Based on similar functions by Jon Yearsley and Hristo Zhivomirov
% Written by Marcin Konowalczyk
% Timmel Group @ Oxford University
 
%% Parse the input
narginchk(0,Inf); nargoutchk(0,1);
 
if nargin < 2 || isempty(beta); beta = 0; end % Default to white noise
assert(isnumeric(beta) amp;amp; isequal(size(beta),[1 1]),'''beta'' must be a number');
assert(-6 <= beta amp;amp; beta <= 6,'''beta'' out of range'); % Put on reasonable bounds
 
%% Generate N-dimensional white noise with 'randn'
X = randn(varargin{:});
if isempty(X); return; end; % Usually happens when size vector contains zeros
 
% Squeeze prevents an error if X has more than one leading singleton dimension
% This is a slight deviation from the pure functionality of 'randn'
X = squeeze(X);
 
% Return if white noise is requested
if beta == 0; return; end;
 
%% Generate corresponding N-dimensional matrix of multipliers
N = size(X);
% Create matrix of multipliers (M) of X in the frequency domain
M = []; 
for j = 1:length(N)
    n = N(j);
    
    if (rem(n,2)~=0) % if n is odd
        % Nyquist frequency bin does not show up in odd-numbered fft
        k = ifftshift(-(n-1)/2:(n-1)/2);
    else
        k = ifftshift(-n/2:n/2-1);
    end
    
    % Spectral multipliers
    m = (k.^2)';
    
    if isempty(M);
        M = m;
    else
        % Create the permutation vector
        M_perm = circshift(1:length(size(M)) 1,[0 1]);
        % Permute a singleton dimension to the beginning of M
        M = permute(M,M_perm);
        % Add m along the first dimension of M
        M = bsxfun(@plus,M,m);
    end
end
% Reverse M to match X (since new dimensions were being added form the left)
M = permute(M,length(size(M)):-1:1);
assert(isequal(size(M),size(X)),'Bad programming error'); % This should never occur
 
% Shape the amplitude multipliers by beta/4 which corresponds to shaping the power by beta
M = M.^(beta/4);
 
% Set the DC component to zero
M(1,1) = 0;
 
%% Multiply X by M in frequency domain
Xstd = std(X(:));
Xmean = mean(X(:));
X = real(ifftn(fftn(X).*M));
 
% Force zero mean unity standard deviation
X = X - mean(X(:));
X = X./std(X(:));
 
% Restore the standard deviation and mean from before the spectral shaping.
% This ensures the random sample from randn is truly random. After all, if
% the mean was always exactly zero it would not be all that random.
X = X   Xmean;
X = X.*Xstd;
end
  

Вот мое решение:

Мой код «min / max» (строки 6 и 7) был плохим. Я хотел разделить все значения в матрице на единственное наибольшее значение в матрице, чтобы все значения были между 0 и 1. Поскольку я неправильно использовал max(), я переходил через максимальное значение каждого столбца и использовал его в качестве своего делителя; таким образом, вертикальные полосы.

В конце концов, вот как выглядит мой код. X — это трехмерная матрица:

 minVal = min(X,[],'all'); % Get the lowest value in the entire matrix
X = X - minVal; % Set min value to zero
maxVal = max(X,[],'all'); % Get the highest value in the entire matrix
X = X ./ maxVal; % Set max value to one