Сетка разбросанных данных без интерполяции

Question

У меня есть три вектора: X, Y и Z, которые представляют трехмерные координаты некоторых данных, найденных в матрице V (т.е. V = f(X,Y,Z)). Ниже приведены некоторые особенности этих данных (предположения / утверждения, если хотите):

• X и Y имеют одинаковую длину, которая (обычно) отличается от Z.

  numel(X) == numel(Y);
  

• Количество уникальных элементов в X обычно отличается от количества уникальных элементов в Y.

• Все значения в Z являются уникальными.

  numel(unique(Z)) == numel(Z);
  

• V имеет размер [numel(X), numel(Z)].

В прошлом я не делал различий между X и Y, и я ссылался на них, используя один индекс, который облегчал построение диаграмм, которые имеют "номер комбинации XY" (то есть 1:numel(X)) на одной оси и Z на другой, чтобы размер V получился удачным.

В настоящее время я хотел бы видеть эффекты X и Y отдельно, поэтому я хочу визуализировать это в 3d, используя смесь scatter3 и contourslice.

Часть рассеяния достаточно проста, так как я могу repmat X и Y вдоль их одноэлементного измерения numel(Z) раз, и аналогично для Z (используя numel(X)). В результате:

Что касается контуров, то для них требуется, чтобы данные были представлены в виде трехмерных массивов, что означает, что я должен размещать точки в структурированной сетке. Создать координаты сетки XX, YY, ZZ просто, используя meshgrid:

[ XX, YY, ZZ ] = meshgrid( unique(X), unique(Y), unique(Z) );

С чем я борюсь, так это с созданием массива 3d VV.

Из рисунка выше мы видим, что большая часть тома не содержит точек - и я бы очень хотел, чтобы оставил это так . Другими словами, идеал VV должен содержать только точки, которые соответствуют исходному набору данных, а остальная часть пространства должна содержать NaN с.

Такие функции, как griddata и interp3 выполняют интерполяцию, которая «заполняет дыры» внутри облака точек - что крайне нежелательно . Я думаю, что индексирование может использоваться здесь для заполнения VV с использованием значений из V, но я не могу придумать, как это сделать.

У меня такой вопрос: как мне сгенерировать VV, который не содержит интерполированных данных ?

Вот минимальный пример:

%% Generate some data:
X = randi(10,100,1);
Y = randi(15,100,1) - 5;
Z = 1:50;
V = X./Y.*Z;

%% Scatter plot:
nXY = numel(X); nZ = numel(Z);
figure();
scatter3( reshape( repmat(X,[1, nZ]),  [], 1), ...
          reshape( repmat(Y,[1, nZ]),  [], 1), ...
          reshape( repmat(Z,[nXY, 1]), [], 1), ...
          [], V(:), '.');

%% Contours:
% Create the 3d grid:
[XX, YY, ZZ] = meshgrid( unique(X), unique(Y), unique(Z) );

% Preallocate VV:
VV = NaN(size(XX));

% Populate VV: <--------------------------------------------- Help needed with this stage
ind = randperm( numel(XX), numel(V) ); % PLACEHOLDER 
VV(ind) = V;

% Plot:
hold on; contourslice(XX, YY, ZZ, VV, X(2), Y(3), Z(10) );

Answer 1

Я полагаю, что возможен излишний путь, если сопоставить все ваши привязанные к сетке индексы со всеми вашими точками линейных данных. Для этого нам нужно ввести несколько измерений, чтобы сравнить 3d-массивы XX и т. Д. С 2d-массивами X и т. Д. Поэлементно:

Xbc = reshape(X, [1,1,1,size(X)]);
Ybc = reshape(Y, [1,1,1,size(Y)]);
Zbc = reshape(Z, [1,1,1,size(Z)]);

Эти массивы изменяются таким образом, что они транслируются с массивами XX и т. Д. Размером [N,M,K] («bc» означает трансляцию). Итак, поэлементное сравнение теперь работает:

match = reshape((XX == Xbc) & (YY == Ybc) & (ZZ == Zbc), [size(XX), numel(V)]);

Этот логический массив имеет размер [N,M,K,P,Q], если V имеет размер [P,Q]. Он содержит ровно столько true с, сколько вы хотите:

>> sum(match(:)) == numel(V)

ans =

  logical

   1

Итак, теперь нам нужно выбрать соответствующие индексы по первым трем измерениям и соединить их с правильным элементом V. Нам нужна смазка для линейных и многоиндексных колен:

[ii,jj,kk,ll] = ind2sub(size(match), find(match));

Теперь все массивы с левой стороны имеют размер [numel(V), 1]; первые три дают вам индексы в XX и т. д., а последний дает индексы в V.

V_inds = ll;
VV_inds = sub2ind(size(VV), ii, jj, kk);

VV(VV_inds) = V(V_inds);

Теперь почему-то я вижу только 3750 уникальных индексов среди 5000 в результате:

>> numel(VV_inds)           

ans =

        5000

>> numel(unique(VV_inds))

ans =

        3750

Я не могу найти другую причину для этого, кроме того, что некоторые из ваших исходных точек данных перекрываются из-за повторений в значениях X и Y, поэтому вы не можете фактически представить их в 3d сетка уникальных точек (потому что некоторые 3d точки содержат более одной точки данных). Я считаю, что следующее доказывает это:

>> size(unique([X,Y], 'rows'))

ans =

    75     2

>> size([X,Y])

ans =

   100     2

Есть 100 (x,y) пар, но только 75 уникальных. Независимо от того, как вы комбинируете их с ортогональными наборами z точек, вы получите повторения в своих точках. Таким образом, вы либо должны отбрасывать избыточность в ваших необработанных данных, либо вам нужно найти другое представление (или взять среднее значение для конфликтующих значений).

---

Я думаю, что у меня есть и более эффективная версия, использующая индексы, сгенерированные unique во время ее запуска. Обратите внимание, что я предполагаю, что вы используете meshgrid вместо ndgrid для генерации сеток, так что размеры результирующих массивов (а также VV) соответствуют уникальным размерам вдоль X, Y и Z соответственно.

% take the indices
[uX, ~, iX] = unique(X);
[uY, ~, iY] = unique(Y);
[uZ, ~, iZ] = unique(Z);

% generate mesh and allocate result
[XX, YY, ZZ] = ndgrid(uX, uY, uZ);
VV = NaN(size(XX));

% switch from `iX`, `iY` and `iZ` to a 2d mesh of size `[P,Q]` where `iX` and `iY` are of size `[P,1]` and `iZ` is of size `[Q,1`]:
% a.k.a. lazy repmat
iXbig = iX + 0*iZ.';
iYbig = iY + 0*iZ.';
iZbig = iZ.' + 0*iX;

% turn 3d indices into linear index into VV
VV_inds = sub2ind(size(VV), iXbig, iYbig, iZbig);

% profit
VV(VV_inds) = V;