CloudRunner

... share your algorithms as a web service

Algorithm: Random Hypersurface Model

Description:
(c) by Marcus Baum & Robin Sandkuehler @ ISAS/KIT
Tags: plot
Usage: Algorithm is public.
Viewed 4557 times, called 532 times
Upload:
Full star Full star Full star Full star Full star
1 vote
Robin Sandkühler
10/13/2012 5:10 p.m. (version #2)

Select Version

Version #2: 10/13/2012 5:10 p.m.
No changelog available

Version #1: 06/05/2012 9:06 p.m.
No changelog available

Run Algorithm

:
: Allow cached result?
Close

Please Wait

Computation is running...

Plots

showing plot # of

Result Value(s)


		

Resource Usage

Execution of the algorithm took seconds.
Peak memory usage for running the algorithm was kilobytes (total usage: kilobytes).

Run ID / Link

ID of this run:
Link:

Result from Cache

Result has been delivered from cache and computed on (UTC).

Matlab log

(show)

		

Error

Using this algorithm in your local MATLAB environment is easy: Click here for instructions!

Usage Instructions for CloudRunner Client

  1. Download the CloudRunner Client by clicking here and place the downloaded file in your MATLAB working directory.

  2. Inside MATLAB, initialize the CloudRunner Client by calling CloudRunner:
    >> CloudRunner

    A login dialog will prompt for your CloudRunner mail address and password. For a start, you can leave the dialog empty and just click "Connect".

    Alternatively, you can provide the login credentials (or empty strings to skip login) as a parameter and hence skip the login dialog. This is useful when using CloudRunner in non-interactive scripts.
    >> CloudRunner('mail@example.com', 'password')

  3. Select this algorithm by its URL. Selecting an algorithm creates the lambda function that proxies calls to the algorithm to the server for execution:
    >> CloudRunnerSelect('http://www.cloudrunner.eu/algorithm/12/random-hypersurface-model/version/2/')

    For the sake of convenience, you can also use the algorithm ID instead of its URL for public algorithms.

  4. Call functions from the algorithm like any regular local function.

Note: You can find further information on the help page.

Source Code

File:

  1 % (c) by Marcus Baum & Robin Sandkuehler @ ISAS/KIT
  2 function randomHypersurfaceModel(numberOfMeasurement)
  3 
  4 % Number of Fourier coefficients
  5 nr_Fourier_coeff = 11;
  6 
  7 % State describtion prior [b0--bn, x, y]
  8 x = zeros(nr_Fourier_coeff + 2, 1);
  9 x(1) = 1.5;
 10 
 11 % State covariance prior
 12 C_x = diag([ones(1, nr_Fourier_coeff).*0.02, 0.3, 0.3]);
 13 
 14 % Measurement noise
 15 measurementNoise = diag([0.2, 0.2].^2);
 16 
 17 % Scale properties
 18 scale.mean = 0.7;
 19 scale.variance = 0.08;
 20 
 21 % Shape resolution for plotting
 22 phi_vec = [0:0.01:2*pi];
 23 
 24 % Object size
 25 a = 3;      % -- width of the horizontal rectangle
 26 b = 0.5;    % | height of the horizontal rectangle
 27 c = 2;      % | height of the vertical rectangle
 28 d = 0.5;    % -- width of the vertical rectangle
 29 
 30 sizeObject = [a b c d];
 31 
 32 % Object shape bounds
 33 objectBounds = [[-d, -c];[d, -c];[d, -b];[a, -b];[a, b];[d, b];[d, c];
 34     [-d, c];[-d, b];[-a, b];[-a, -b];[-d, -b]]' ./ 2;
 35 
 36 
 37 %% Main
 38 
 39 % Plot
 40 h_object = fill(objectBounds(1, :), objectBounds(2, :), [.7 .7 .7]);
 41 hold on
 42 xlim([-3 3]);
 43 ylim([-3 3]);
 44 axis equal
 45 xlabel('x-Axis')
 46 ylabel('y-Axis')
 47 title('Random Hypersurface Model Simulation')
 48 
 49 
 50 
 51 for j = 1 : numberOfMeasurement
 52     
 53     % Get new measurement
 54     newMeasurement = getNewMeasurement(sizeObject, measurementNoise);
 55     
 56     % Filter step
 57     [x, C_x] = UKF_FilterStep(x, C_x, newMeasurement, [scale.mean; [0 0]'], ...
 58         blkdiag(scale.variance, measurementNoise), @f_meas_pseudo_squared, nr_Fourier_coeff);
 59     
 60     
 61     % Plot
 62     shape = calcShape(phi_vec, x, nr_Fourier_coeff);
 63     
 64     h_measure = plot(newMeasurement(1), newMeasurement(2), '+');
 65     h_shape =  plot(shape(1, :), shape(2, :), 'g-', 'linewidth', 2);
 66     legend([h_object, h_measure, h_shape],'Target', 'Measurement', 'Estimated shape')
 67     drawnow;
 68     
 69     if j ~= numberOfMeasurement
 70         delete(h_shape)
 71     end
 72     
 73 end
 74 
 75 hold off
 76 end
 77 
 78 
 79 
 80 function measurement =  getNewMeasurement( sizeObject, measurementNoise )
 81 
 82 a = sizeObject(1); % -- width of the horizontal rectangle
 83 b = sizeObject(2); % | height of the horizontal rectangle
 84 c = sizeObject(3); % | height of the vertical rectangle
 85 d = sizeObject(4); % -- width of the vertical rectangle
 86 
 87 measurementsourceNotValid = 1;
 88 
 89 while measurementsourceNotValid
 90     
 91     %Measurementsoure
 92     x = -a/2 + a.*rand(1, 1);
 93     y = -c/2 + c.*rand(1, 1);
 94     
 95     if y > b/2 && x < -d/2 || y > b/2 && x > d/2 ||...
 96        y < -b/2 && x < -d/2 || y < -b/2 && x > d/2
 97         
 98         x = -a/2 + a.*rand(1, 1);
 99         y = -c/2 + c.*rand(1, 1);
100         
101     else
102         measurementsourceNotValid = 0;
103     end
104 end
105 
106 % Add zero-mean gaussian noise to the measurement sources
107 measurement = [x; y] + (randn(1, 2) * chol(measurementNoise))';
108 end
109 
110 function pseudoMeasurement = f_meas_pseudo_squared(x, noise, y, nr_Fourier_coeff)
111 
112 numberOfSigmaPoints = size(x, 2);
113 pseudoMeasurement = zeros(1, numberOfSigmaPoints);
114 
115 for j = 1 : numberOfSigmaPoints
116     
117     s = noise(1, j);
118     v = noise(2:3, j);
119     b = x(1:nr_Fourier_coeff, j);
120     m = x(nr_Fourier_coeff + 1:nr_Fourier_coeff + 2, j);
121     
122     theta = atan2(y(2) -m(2), y(1) - m(1))+2*pi;
123 
124     
125     R = calcFourierCoeff(theta, nr_Fourier_coeff);
126     
127     e = [cos(theta); sin(theta)];
128     
129     pseudoMeasurement(j) = (norm( m - y ))^2 - (s^2 *(R * b).^2 + 2 * s * R * b * e' * v + norm(v)^2);
130 end
131 end
132 
133 
134 
135 
136 function fourie_coff = calcFourierCoeff(theta, nr_Fourier_coeff)
137 
138 fourie_coff(1) = 0.5;
139 
140 index = 1;
141 for i = 2 : 2 : (nr_Fourier_coeff - 1)
142     fourie_coff(i : i + 1) = [cos(index * theta) sin(index * theta)];
143     index = index + 1;
144 end
145 
146 
147 end
148 
149 
150 function shape = calcShape(phi_vec, x, nr_Fourier_coeff)
151 
152 shape = zeros(2, length(phi_vec));
153 
154 for i = 1 : length(phi_vec)
155     phi = phi_vec(i);
156     R = calcFourierCoeff(phi, nr_Fourier_coeff);
157     e = [cos(phi) sin(phi)]';
158     shape(:, i) = R * x(1:end - 2) * e + x(end - 1:end);
159 end
160 
161 end
162 
163 function [x_e, C_e] = UKF_FilterStep(x, C, measurement, measurementNoiseMean, measurementNoiseCovariance, measurementFunctionHandle, numberOfFourierCoef)
164 
165 alpha = 1;
166 beta = 0;
167 kappa = 0;
168 
169 %% Calculate Sigma Points
170 
171 %Merge State and noise mean
172 x_ukf = [x; measurementNoiseMean];
173 
174 
175 %Merge State and noise Covariance
176 C_ukf = blkdiag(C, measurementNoiseCovariance);
177 
178 n = size(x_ukf, 1);
179 n_state = size(x, 1);
180 
181 lamda = alpha^2 * (n + kappa) - n;
182 
183 
184 % Calculate Weights Mean
185 WM(1) = lamda / (n + lamda);
186 WM(2 : 2 * n + 1) = 1 / (2 * (n + lamda));
187 % Calculate Weights Covariance
188 WC(1) = (lamda / (n + lamda)) + (1 - alpha^2 + beta);
189 WC(2 : 2 * n + 1) = 1 / (2 * (n + lamda));
190 
191 %Calculate Sigma Points
192 A = sqrt(n + lamda) * chol(C_ukf)';
193 
194 xSigma = [zeros(size(x_ukf)) -A A];
195 xSigma = xSigma + repmat(x_ukf, 1, size(xSigma, 2));
196 
197 %% Filterstep
198 z = 0;
199 C_yy = 0;
200 C_xy = 0;
201 
202 zSigmaPredict = feval(measurementFunctionHandle, xSigma(1:n_state,:), xSigma(n_state + 1:n, :), measurement, numberOfFourierCoef );
203 
204 for i = 1 : size(zSigmaPredict, 2);
205     z = z + ( zSigmaPredict(:,i) * WM(i) );
206 end
207 
208 
209 for i = 1 : size(zSigmaPredict, 2)
210     C_yy = C_yy + WC(i) * ( (zSigmaPredict(:,i) - z ) * ( zSigmaPredict(:,i) - z )') ;
211     C_xy = C_xy + WC(i) * ( (xSigma(1:size(x, 1),i) - x ) * ( zSigmaPredict(:,i) - z )');
212 end
213 
214 K = C_xy / C_yy;
215 x_e = x + K * (zeros(size(z)) - z);
216 C_e = C - K * (C_yy) * K';
217 
218 
219 end
Download algorithm (1 file) as ZIP

Comments

Please login to post a comment.

#2: Chaoqun Yang on 03/10/2016 4:23 a.m.
wonderful

#1: Christian Mandery on 06/11/2012 6:49 p.m.
Awesome plots!