smooth scattered data of a multidimensional matrix

WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

Nesta

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 30.07.2013, 19:05 Titel: smooth scattered data of a multidimensional matrix

Hallo zusammen,

ich habe eine multidimensionale Matrix mit gescatterten Daten, die ich fitten/smoothen möchte. Für zwei Dimensionen habe ich bereits eine Lösung:

x und y sind Vektoren mit den Koordinaten meiner Daten.
X und Y sind dann Matrizen mit einem Koordinatengitter.
In Z habe ich bereits die gescatterten Daten über diesem zwei dimensionalen Gitter. Die Funktion gridfit smootht mir dann die verrauschten Daten aus Z in die Matrix F. x_len und y_len geben die Knoten der Interpolationspunkte an.

Code:

[X,Y] = ndgrid(x , y);

[F,X,Y] = gridfit( X, Y, Z, x_len, y_len, 'smoothness', 15);

Funktion ohne Link?

X, Y, Z und F sind also 2D-Matizen gleicher Größe! Die stark zackigen Contourlinen, die mir contour(X, Y, Z) lieferen würde, werden mir contour(X, Y, F) so beliebig gesmootht.

gridfit ist aber leider nur für 2D-Matizen. Ich suche eine Verallgemeinerung, also eine Funktion oder Methode, die mir etwas ähnlihches für 1 oder >= 3 Dimensionen macht.

Hat jemand vielleicht einen Hinweis ob es da schon eine fertige Funktion gibt, oder ein Schlagwort nach dem ich suchen könnte.

Vielen Dank schon mal im voraus.

Nesta

Thomas84

Forum-Meister


	Beiträge: 546

	Anmeldedatum: 10.02.10

	Wohnort: ---

	Version: ---

Verfasst am: 31.07.2013, 08:37 Titel:

Hallo Nesta,

ich hatte mal das gleiche Problem und hab dafür mal eine Funktion geschrieben. Vielleicht hilft sie dir bei deiner Aufgabe:

Code:

function [V_grid,hess,varargout] = scattered2grid(Xi,V,nodes,method,smooth)
% scattered2grid
%
% Ermittlung eines nd-Kennfeldes dasierend auf unregelmäßig verteilten Daten.
% Dabei wird grob gesprochen die Norm ||Xi - interpne(V,Xi,nodelist)|| + ||d^2 V_grid/d Xi^2|| minimiert.
%
% Der erste Term ist der Abstand zwischen den Messdaten und den aus dem
% Kennfeld interpolierten Werten. Er beschreibt demnach wie gut das
% Kennfeld und die Messwerte überseinstimmen.
% Der zweite Term beschreibt die zweite Ableitung und führt
% zu einer Glättung des Kennfeldes.
%
% Syntax:
% [V_grid,hess,X1,...,Xn] = scattered2grid(Xi,V,nodelist,method,smooth)
%
% Argumente:
% Xi - n x p array mit unregelmäßig verteilten Daten
% V - n x 1 array mit Funktionswerten
% nodelist - 1 x p cell array mit Stützstellen für jede Dimension
% method - string mit folgenden möglichen Inhalt: nearest {linear} spline pchip cubic
% (Beschreibung siehe interp1)
% smooth - p x 1 array mit Parameter der angibt wie stark das Kennfeld in Richtung der
% jeweilligen Dimension geglättet werden soll (große Werte
% führen zu starker Glättung).
% Falls nur ein Wert angegeben wird, wird der gleiche Wert für alle Dimensionen angenommen.
%
% Rückgabewerte:
% V_grid - p dimensionales array mit Kennfeld (Dimensionen entsprechen der Länge der Vektoren in nodelist)
% hess - Hessematrix (Kovarianz der ermittelten Parameter = 0.5*inv(hess)). Elemente mit Werten
% nahe Null können nicht ermittelt werden
% X1,...,Xn - Stützstellen des Kennfeld (generiert durch nd-grid)
%
% Beispiel 2d:
% x = rand(100,1);
% y = rand(100,1);
% z = exp(x+2*y);
% nodes = {0:.1:1,0:.1:1};
%
% g = gridfit(x,y,z,nodes{1},nodes{2});
% [s,hess,xx,yy] = scattered2grid([x,y],z,nodes,'linear',0.01);
%
% figure('name','Kennfeld');
% % y und x müssen vertauscht werden, da sie durch nd-grid und nicht durch
% % meshgrid erstellt werden
% stem3(yy,xx,g);hold on;
% stem3(xx,yy,s,'r');
% xlabel('x');ylabel('y');zlabel('z = e^{x+2y}');
% legend('gridfit','scattered2grid');
%
% figure('name','Interpolation');
% stem3(x,y,z);hold on;
% stem3(x,y,interpne(s,[x,y],nodes),'g');
% stem3(x,y,interpne(g.',[x,y],nodes),'r');
% legend('original','scatterde2grid','gridfit');
%
% Beispiel 3d:
% % Kennfeld mit 216 Stützstellen wird anhand von nur 100 Messwerten berechnet
% nodes = {0:0.2:1,0:0.2:1,0:0.2:1};
% Xi = rand(100,3);
% V = tanh(Xi*[1.2;2.3;-3.4]);
% [V_grid,h,x,y,z] = scattered2grid(Xi,V,nodes,'linear',0.01);
%
% figure('name','scattered2grid - 3d');rotate3d on;
% for k = 1:length(z)
% subplot(3,2,k);
% surf(x(:,:,k),y(:,:,k),V_grid(:,:,k));hold on;
% stem3(x(:,:,k),y(:,:,k),tanh(1.2*x(:,:,k) + 2.3*y(:,:,k) - 3.4*nodes{3}(k)));
% xlabel('x');ylabel('y');title(['z = ',num2str(nodes{3}(k))]);
% end
%
% see also: gridfit, interp1, interpne, ndgrid
%
% ------
% Author: Thomas Lehmann
% Created: 2012-12-04

if nargin < 4
method = 'linear';
end

if nargin < 5
smooth = 1;
end

n = size(Xi,1); % Anzahl der Datenpunkte
dim = size(Xi,2); % Anzahl der Dimensionen (== length(nodes)
sz = cellfun(@length,nodes); % Size
s = prod(sz); % Anzahl der Stützpunkte

if length(nodes) ~= dim
error('Anzahl der Stützstellenvektoren muss mit Dimension übereinstimmen');
end

if numel(smooth) ~= dim
smooth = smooth*ones(dim,1);
end

%% interne Skalierung der Stützstellenvektoren
org_nodes = nodes;
for k = 1:length(nodes)
nodes{k} = (nodes{k}-nodes{k}(1))/(nodes{k}(end)-nodes{k}(1));
Xi(:,k) = (Xi(:,k)-org_nodes{k}(1))/(org_nodes{k}(end)-org_nodes{k}(1));
end

%% A Matrix

for k = 1:dim
pp{k} = interp1(nodes{k},eye(length(nodes{k})),method,'pp');
phi{k} = ppval(pp{k},Xi(:,k));

phi{k}(phi{k} < 0) = 0;
phi{k}(phi{k} > 1) = 1;
end

A = zeros(n,s);
temp = zeros(n,dim);

for k = 1:s
ind = myind2sub(sz,k);
for j = 1:dim
temp(:,j) = phi{j}(:,ind(j));
end
A(:,k) = prod(temp,2);
end

%% Anzahl Bedingungen (2.Ableitung = 0)
n_m = zeros(1,dim);
for l = 1:dim
temp = sz;
if temp(l) > 2
temp(l) = temp(l) - 2;
n_m(l) = prod(temp);
else
n_m(l) = 0;
end
end

n_bed = sum(n_m);

%% Gewichtung berechnen

for j = 1:dim
deltap = [diff(nodes{j}),NaN];
deltam = [NaN,-diff(nodes{j})];
t1 = 2./(deltap.^2+deltam.^2);
t2 = (deltap + deltam)./(deltap - deltam);

gm{j} = t1 + t1.*t2;
g0{j} = -2*t1;
gp{j} = t1 - t1.*t2;
end

%% B berechnen
B = zeros(n_bed,s);

j = 1;
for k = 1:s
sub_k = myind2sub(sz,k);
for l = 1:dim
if sub_k(l) > 1 && sub_k(l) < sz(l)
v = zeros(1,dim);
B(j,mysub2ind(sz,sub_k + v)) = g0{l}(sub_k(l))*smooth(l);
v(l) = 1;
B(j,mysub2ind(sz,sub_k + v)) = gp{l}(sub_k(l))*smooth(l);
v(l) = -1;
B(j,mysub2ind(sz,sub_k + v)) = gm{l}(sub_k(l))*smooth(l);
j = j + 1;
end
end
end

V_grid = [A;B]\[V;zeros(n_bed,1)];

V_grid = reshape(V_grid,cellfun(@length,nodes));

% Hesse-Matrix verwenden um zu ermitteln welche Stützstellen ermittelt
% werden konnten. Die restlichen auf NaN setzen.
hess = 2*[A;B].'*[A;B];
ind = find(diag(hess) == 0);
V_grid(ind) = NaN;

%% ndgrid
nodes = org_nodes;
if dim==1, nodes = repmat(nodes,[1 max(nargout,2)]); end
nin = length(nodes);
nout = max(nargout-2,nin);

for i=length(nodes):-1:1,
nodes{i} = full(nodes{i}); % Make sure everything is full
siz(i) = numel(nodes{i});
end
if length(siz)<nout, siz = [siz ones(1,nout-length(siz))]; end

varargout = cell(1,nout);
for i=1:nout,
x = nodes{i}(:); % Extract and reshape as a vector.
sn = siz; sn(i) = []; % Remove i-th dimension
x = reshape(x(:,ones(1,prod(sn))),[length(x) sn]); % Expand x
varargout{i} = permute(x,[2:i 1 i+1:nout]); % Permute to i'th dimension
end

function vout = myind2sub(siz,ndx)
% sub = myind2sub(siz,lin_ind)
% wie ind2sub nur das array zurückgegeben wird anstatt einzelne skalare
%
% Argumente:
% siz - 1 x dim Array mit Dimension des Arrays
% lin_ind - n x 1 Array mit linearen Indizes
%
% Rückgabewerte:
% sub - n x dim Array mit subscript indizes
%
% Beispiel:
% sz = [6,4,9];
% ind = [34,44,47,48,73];
% [I1,I2,I3] = ind2sub(sz,ind)
% sub = myind2sub(sz,ind)
%
% see also: ind2sub, mysub2ind

if numel(ndx) > 1
vout = zeros(numel(ndx),length(siz));
for j = 1:numel(ndx)
vout(j,:) = myind2sub(siz,ndx(j));
end
return
elseif numel(ndx) == 0
vout = [];
return
end

nout = length(siz);
siz = double(siz);

vout = zeros(size(siz));

if length(siz)<=nout,
siz = [siz ones(1,nout-length(siz))];
else
siz = [siz(1:nout-1) prod(siz(nout:end))];
end
n = length(siz);
k = [1 cumprod(siz(1:end-1))];
for i = n:-1:1,
vi = rem(ndx-1, k(i)) + 1;
vj = (ndx - vi)/k(i) + 1;
vout(i) = vj;
ndx = vi;
end

function ndx = mysub2ind(siz,sub)
% lin_ind = mysub2ind(sz,sub)
%
% wie sub2ind nur mit anderem Format des Argumentes sub
%
% Argumente:
% size - 1 x dim Array mit Dimension des Arrays
% sub - n x dim Array mit subscript indizes
%
% Rückgabewerte:
% lin_ind - n x 1 Array mit linearen Indizes
%
% Beispiel:
% sz = [6,4,9];
% ind = [34,44,47,48,73];
% [I1,I2,I3] = ind2sub(sz,ind);
% sub = myind2sub(sz,ind);
% ind = mysub2ind(sz,sub)
%
% see also: sub2ind, myind2sub

if isempty(sub)
ndx = [];
return;
end

siz = double(siz);

dim = size(sub,2);
len = size(sub,1);

if length(siz) ~= dim
error('MATLAB:sub2ind:InvalidSize','Length of siz must match size(sub,2)');
end

%Compute linear indices
k = [1 cumprod(siz(1:end-1))];
ndx = ones(len,1);

for l = 1:len
for i = 1:dim
v = sub(l,i);
if v < 1 || v > siz(i)
%Verify subscripts are within range
ndx(l) = NaN;
break
end
ndx(l) = ndx(l) + (v-1)*k(i);
end
end

Funktion ohne Link?

viele Grüße
Thomas

Nesta

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 31.07.2013, 14:07 Titel:

Hallo Thomas,

vielen Dank für deine Hilfe. Ich hab deinen Code mal ausprobiert und im Prinzip ist das auch was ich suche. Das Problem ist nur, dass ich in jede Diminsion so zwischen 100 und 200 Interpolationspunkte habe, und dann wird die Hessematrix einfach zu groß. Selbst in nur 2 Dimensionen hat es zu lange gedauert, bzw. war dann sogar irgendwann der Speicher voll.

Ich hab inzwischen zwei Funktionen gefunden, die genau das machen was ich brauche (ohnen Hessematrix). Beide heißen "smoothn" und sind zu finden unter:
http://www.mathworks.com/matlabcentral/fileexchange/

meine Empfehlung:
http://www.mathworks.com/matlabcent.....thing-for-1-d-to-n-d-data

Trotzdem nochmal vielen Dank

Nesta


Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:	Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen

Impressum | Nutzungsbedingungen | Datenschutz | FAQ |

RSS

Hosted by:

Copyright © 2007 - 2025 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.