Abbreviations in JAG3D's FormFittingToolbox

Hi, I used your nice software ( JAG3D.exe) for fitting cone and others object. From the statistical report after fitting object (e.g cone, sphere and plane)
Can I ask what rx, ry, rz, and how T-post and T priori have been calculated for each point after fitting. Also are these (σx, σy, and σz)mean standard deviation of a number of iteration for each point.
I would be happy if you could explain each elements in the statistical result of fitting object such as (rx,ry,rz; T-priori, T-post; σx, σy, and σz )

Regards

Eintrag gesperrt

16273 Views

Hi Balata,

Can I ask what rx, ry, rz,

r{x,y,z} denotes the redundancy of the coordinate component of each point. The sum is equal to the degree of freedom of the least squares problem. The value is defined between $0 \le r \le 1$ (or similar 0-100%).

and how T-post and T priori have been calculated for each point after fitting.

Both values are derived by a multiple test. It is a generalization of Baarda's data-snooping. A description is published by e.g. Jäger et al. or Koch. The formulas are also given in equation (16) and (17).

Also are these (σx, σy, and σz) mean standard deviation of a number of iteration for each point.

I'm not sure what do you mean with number of iteration for each point. With σ{x,y,z} the uncertainties of a point is given. These values are taken from the estimated variance-covariance-matrix of the observations.

Kind regards
Micha

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Tags:
JAG3D, FormFittingToolbox, Abbreviation

Eintrag gesperrt

16241 Views

Hi,

Sorry for being late on replying to your nice solution, It is my fault I though that I will receive your nice answer through the web. Thank you very much for responding my question and for your effort.

could you a little bit more explain about rx,ry and rx please? Do you mean the sum of redundancy of the coordinate components of each point is equal to the degree of freedom of the least squares problem ( is the number of iteration used for computing each component until convergence achieved, and then the sum of is equal to the degree of freedom?). or for e.g if I have 182 pts captured by the laser scanner on the surface of cone. The degree of freedom will be 182 -7parameter =175 redundancy.. . I little bit confused)

Regarding the σ{x,y,z} I meant what are represent in the statistic table for each point.

Sorry to bother you, I have another question which is so important is the pyramid object. Can the software JAG3D able to fit the pyramid and calculate the statistical result of its apex as it does for cone, sphere and others. if not what is the solution

I have a set of 3D data captured by Laser Scanner on a pyramid 4faces. I attached a copy of points here.if you have any question please let me know.
Again many thanks for you, I really thanks a lot your contribute and your answer

-8.774 10.2929 -0.2645
-8.7546 10.2701 -0.3177
-8.8042 10.2864 -0.2647
-8.7602 10.2767 -0.2908
-8.7154 10.3041 -0.2639
-8.7347 10.2841 -0.2905
-8.747 10.2985 -0.2641
-8.8216 10.3067 -0.2113
-8.8498 10.2977 -0.2114
-8.7888 10.3102 -0.2109
-8.8108 10.294 -0.2379
-8.7844 10.3051 -0.2378
-8.7524 10.305 -0.2374
-8.7637 10.3184 -0.2109
-8.8621 10.312 -0.1847
-8.828 10.3142 -0.1843
-8.7998 10.3231 -0.1842
-8.7734 10.3298 -0.1841
-8.7055 10.3322 -0.2103
-8.7251 10.3155 -0.2372
-8.6989 10.3245 -0.2371
-8.7321 10.3238 -0.2106
-8.6751 10.3405 -0.2103
-8.7436 10.3375 -0.1838
-8.7142 10.3427 -0.1837
-8.6877 10.3555 -0.1837
-8.6584 10.3609 -0.1835
-8.9011 10.314 -0.4292
-8.8979 10.3103 -0.4021
-8.8503 10.2982 -0.4009
-8.8934 10.3051 -0.3204
-8.8922 10.3037 -0.3744
-8.8957 10.3077 -0.3475
-8.8903 10.3014 -0.266
-8.8403 10.2866 -0.3192
-8.8439 10.2908 -0.3733
-8.8447 10.2917 -0.3463
-8.7987 10.2798 -0.3454
-8.7972 10.2781 -0.3185
-8.842 10.2886 -0.2651
-8.8414 10.2879 -0.292
-8.7972 10.2781 -0.2914
-8.8914 10.3028 -0.2118
-8.8913 10.3026 -0.239
-8.8911 10.3024 -0.1848
-8.8431 10.2899 -0.2381
-8.7668 10.3218 -0.4267
-8.8744 10.3263 -0.456
-8.7807 10.3382 -0.4543
-8.8118 10.3373 -0.4551
-8.8147 10.2984 -0.4003
-8.8596 10.3091 -0.4282
-8.8263 10.3121 -0.4275
-8.7545 10.3073 -0.3992
-8.7821 10.3024 -0.3996
-8.7932 10.3154 -0.4269
-8.7576 10.2737 -0.3445
-8.8032 10.285 -0.3726
-8.7698 10.288 -0.3719
-8.6894 10.3576 -0.4529
-8.7478 10.3424 -0.4536
-8.7209 10.3506 -0.4534
-8.7354 10.3278 -0.4262
-8.7101 10.3377 -0.4259
-8.7437 10.2946 -0.3717
-8.6503 10.351 -0.3985
-8.6315 10.3681 -0.3985
-8.6334 10.3705 -0.452
-8.6954 10.3203 -0.3985
-8.6694 10.3337 -0.3982
-8.6744 10.3396 -0.4252
-8.655 10.3567 -0.4252
-8.6332 10.3702 -0.4251
-8.7059 10.2929 -0.3174
-8.7298 10.2782 -0.3441
-8.7117 10.2996 -0.3713
-8.7069 10.294 -0.3442
-8.6881 10.3115 -0.3442
-8.7283 10.2765 -0.3173
-8.6858 10.3089 -0.3174
-8.6907 10.3147 -0.3712
-8.6435 10.343 -0.3173
-8.6684 10.3324 -0.3711
-8.6662 10.3298 -0.344
-8.6463 10.3463 -0.3711
-8.6288 10.365 -0.3712
-8.6454 10.3453 -0.3443
-8.6273 10.3631 -0.3443
-8.7063 10.2933 -0.2904
-8.6863 10.3095 -0.2636
-8.6875 10.3109 -0.2905
-8.6632 10.3263 -0.2636
-8.6636 10.3268 -0.2903
-8.6423 10.3415 -0.2903
-8.6206 10.355 -0.2634
-8.622 10.3567 -0.3174
-8.6235 10.3586 -0.2903
-8.6664 10.3301 -0.2369
-8.6401 10.3389 -0.2366
-8.6232 10.3582 -0.2099
-8.6199 10.3542 -0.2366

Eintrag gesperrt

16262 Views

Hi Balata,

could you a little bit more explain about rx,ry and rx please?

It is quite easy. The overall redundancy of your least squares problem - also known as degree of freedom - is the number of unnecessary (reading as redundant) observations. An example: A circle in 2D has u=3 parameters, the radius and the center point represented by the values x0 and y0. To estimate a circle, you need at least n=3 points (which do not lie on a straight line). The degree of freedom is zero: f = n-u = 3-3 = 0. If f is zero, you don't need a least square solver because the derived solution is unique and the residuals v are zero. Thus, a closed solution exists.

A least squares solution is needed, if f > 0 to find the most probable solution with the condition $v^Tv == Min$ . Whereas u, the number of unknown parameters, is constant, the number of observations (or points) n has to increase to enlarge f.

A circle fit, which was carried out with 4 points, has a degree of freedom of f = n-u = 4-3 = 1. The integer f can be broken down to the influence of each point - the redundancy rx, ry, rz. Whereas a value of r=0 means, the observation completely affected the final solution, denotes a value r=1 an unneeded observation. Therefore, r is an important parameter because it gives an impression of the controlled unit. If r=0 it is impossible to detect blunders, on the other hand r=1 means a full controlled point. A recommended value is 0.4 < r < 0.7 - a balance between cost and benefit.

Do you mean the sum of redundancy of the coordinate components of each point is equal to the degree of freedom of the least squares problem

Yes, because trace(R) = n-u, where R is the so-called redundancy matrix. The values rx, ry and rz are taken from the main diagonal of R.

Regarding the σ{x,y,z} I meant what are represent in the statistic table for each point.

The σ{x,y,z} denote the uncertainties _after_ the adjustment. Because of the imperfection of your instrument, each observation gets a (small) residual during the adjustment. If you have a perfect surface (or a perfect mathematical model), the estimated uncertainties demonstrate the accuracy of your instrument.

Can the software JAG3D able to fit the pyramid and calculate the statistical result of its apex as it does for cone, sphere and others.

A pyramid is currently not supported by the FormFittingToolbox because a pyramid is a composite figure of k planes (with k>2). To get the pyramid-apex, you have to estimate the k planes of each side. Finally, the intersection point of the planes is the apex.

If you separate your cloud for each plane, I try to develop a solution in Octave.

Kind regards
Micha

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Eintrag gesperrt

16192 Views

Hi, Thank you very much for your valuable help and support. I have got the point. i was busy with your nice articular (reference point determination with a new mathematical model...).

Actually, I did a program using Matlab code contain two parts: The first Part is automatic detection of the target from surrounding objects in the laser point cloud. The second part is computing the geometric centre of the detected target. for example for the pyramid faces i used Plane intersection algorithm after each face separately detected. The problem is I don't know how to calculate statistically the precision of the apex as your JAG3D did very well with cone and sphere...

How to calculate σxyz of the pyramid apex. What I did so far is σxyz (apex)= σx(face1)+σx(face2)+σx(face3)+σx(face3). This is diffidently incorrect. therefore I would like to help me, how to propagate these error to determine the σapex. I also tried to solve it using inverse plan equation such as [ σa1σx+σb1σy+σc1σz=1; σa2σx+σb2σy+σc2σz=1; σa3σx+σb3σy+σc3σz=1; σa4σx+σb4σy+σc4σz=1]. The 4 set of parameters ( σa1,σb1,σc1 to σa4,σb4,σc4) are known from fitting each face separately using JAG3D and then, what is left is the σx,σy,σz. This method is also I'm not happy with it.

A copy of separate data of each face of the pyramid are as follows

face1
-8.7740 10.2929 -0.2645
-8.7546 10.2701 -0.3177
-8.8042 10.2864 -0.2647
-8.7602 10.2767 -0.2908
-8.7154 10.3041 -0.2639
-8.7347 10.2841 -0.2905
-8.7470 10.2985 -0.2641
-8.8216 10.3067 -0.2113
-8.8498 10.2977 -0.2114
-8.7888 10.3102 -0.2109
-8.8108 10.2940 -0.2379
-8.7844 10.3051 -0.2378
-8.7524 10.3050 -0.2374
-8.7637 10.3184 -0.2109
-8.8621 10.3120 -0.1847
-8.8280 10.3142 -0.1843
-8.7998 10.3231 -0.1842
-8.7734 10.3298 -0.1841
-8.7055 10.3322 -0.2103
-8.7251 10.3155 -0.2372
-8.6989 10.3245 -0.2371
-8.7321 10.3238 -0.2106
-8.6751 10.3405 -0.2103
-8.7436 10.3375 -0.1838
-8.7142 10.3427 -0.1837
-8.6877 10.3555 -0.1837
-8.6584 10.3609 -0.1835
Face2
-8.9011 10.3140 -0.4292
-8.8979 10.3103 -0.4021
-8.8503 10.2982 -0.4009
-8.8934 10.3051 -0.3204
-8.8922 10.3037 -0.3744
-8.8957 10.3077 -0.3475
-8.8903 10.3014 -0.2660
-8.8403 10.2866 -0.3192
-8.8439 10.2908 -0.3733
-8.8447 10.2917 -0.3463
-8.7987 10.2798 -0.3454
-8.7972 10.2781 -0.3185
-8.8420 10.2886 -0.2651
-8.8414 10.2879 -0.2920
-8.7972 10.2781 -0.2914
-8.8914 10.3028 -0.2118
-8.8913 10.3026 -0.2390
-8.8911 10.3024 -0.1848
-8.8431 10.2899 -0.2381
Face 3
-8.7668 10.3218 -0.4267
-8.8744 10.3263 -0.4560
-8.7807 10.3382 -0.4543
-8.8118 10.3373 -0.4551
-8.8147 10.2984 -0.4003
-8.8596 10.3091 -0.4282
-8.8263 10.3121 -0.4275
-8.7545 10.3073 -0.3992
-8.7821 10.3024 -0.3996
-8.7932 10.3154 -0.4269
-8.7576 10.2737 -0.3445
-8.8032 10.2850 -0.3726
-8.7698 10.2880 -0.3719
-8.6894 10.3576 -0.4529
-8.7478 10.3424 -0.4536
-8.7209 10.3506 -0.4534
-8.7354 10.3278 -0.4262
-8.7101 10.3377 -0.4259
-8.7437 10.2946 -0.3717
Face 4
-8.6503 10.3510 -0.3985
-8.6315 10.3681 -0.3985
-8.6334 10.3705 -0.4520
-8.6954 10.3203 -0.3985
-8.6694 10.3337 -0.3982
-8.6744 10.3396 -0.4252
-8.6550 10.3567 -0.4252
-8.6332 10.3702 -0.4251
-8.7059 10.2929 -0.3174
-8.7298 10.2782 -0.3441
-8.7117 10.2996 -0.3713
-8.7069 10.2940 -0.3442
-8.6881 10.3115 -0.3442
-8.7283 10.2765 -0.3173
-8.6858 10.3089 -0.3174
-8.6907 10.3147 -0.3712
-8.6435 10.3430 -0.3173
-8.6684 10.3324 -0.3711
-8.6662 10.3298 -0.3440
-8.6463 10.3463 -0.3711
-8.6288 10.3650 -0.3712
-8.6454 10.3453 -0.3443
-8.6273 10.3631 -0.3443
-8.7063 10.2933 -0.2904
-8.6863 10.3095 -0.2636
-8.6875 10.3109 -0.2905
-8.6632 10.3263 -0.2636
-8.6636 10.3268 -0.2903
-8.6423 10.3415 -0.2903
-8.6206 10.3550 -0.2634
-8.6220 10.3567 -0.3174
-8.6235 10.3586 -0.2903
-8.6664 10.3301 -0.2369
-8.6401 10.3389 -0.2366
-8.6232 10.3582 -0.2099
-8.6199 10.3542 -0.2366

Eintrag gesperrt

16193 Views

Hi,

How to calculate σxyz of the pyramid apex.

I used a combined algorithm which estimates the 4 plane parameters and the apex in one pour, formulated as a Gauß-Helmert model with unknowns in the restrictions. The Octave (quick and dirty) code is shown below. The first three parameters represent the apex and its uncertainties.

         -8.74767399063614       0.00221451170576236
          10.2640449004265        0.0011073336551347
        -0.328875210651509       0.00218131552759802
        -0.210956422930505       0.00890919575377389
         0.873711747603854       0.00478721066144033
        -0.438320852485237       0.00927348301769577
          10.9573474841569        0.0779700931038434
         0.275432828201474        0.0122221368827549
         0.960863952854683       0.00351702329905086
        0.0296179211527723       0.00860389157635439
          7.44321356781815         0.144466915923314
        -0.151048445180385        0.0117018439751209
         0.828459948964864       0.00695058710419911
         0.539294428090744        0.0104652274965152
          9.64731210096934          0.11476352555406
        -0.625323330824657       0.00737061968100696
         0.775035753685848       0.00594490104169702
         0.091051152848087       0.00677300352475538
          13.3951819449872        0.0109245236511967

kind regards
Micha

 
%% 4-Faces-Pyramid
%
%   nx(i) * x(j,i) + ny(i) * y(j,i) + nz(i) * z(j,i) - d(i) = 0
% with i = 1 .. 4
%
%   ||n(i)|| == 1
% 
%  and 
%
%   n(i)' * Apex -d(i) == 0 
%
%
% @version 1.01
% @author Michael Loesler - derletztekick.com
 
clear all
close all
format long g
clc
 
plane1 = [
 -8.7740 10.2929 -0.2645
 -8.7546 10.2701 -0.3177
 -8.8042 10.2864 -0.2647
 -8.7602 10.2767 -0.2908
 -8.7154 10.3041 -0.2639
 -8.7347 10.2841 -0.2905
 -8.7470 10.2985 -0.2641
 -8.8216 10.3067 -0.2113
 -8.8498 10.2977 -0.2114
 -8.7888 10.3102 -0.2109
 -8.8108 10.2940 -0.2379
 -8.7844 10.3051 -0.2378
 -8.7524 10.3050 -0.2374
 -8.7637 10.3184 -0.2109
 -8.8621 10.3120 -0.1847
 -8.8280 10.3142 -0.1843
 -8.7998 10.3231 -0.1842
 -8.7734 10.3298 -0.1841
 -8.7055 10.3322 -0.2103
 -8.7251 10.3155 -0.2372
 -8.6989 10.3245 -0.2371
 -8.7321 10.3238 -0.2106
 -8.6751 10.3405 -0.2103
 -8.7436 10.3375 -0.1838
 -8.7142 10.3427 -0.1837
 -8.6877 10.3555 -0.1837
 -8.6584 10.3609 -0.1835
];
 
plane2 = [
 -8.9011 10.3140 -0.4292
 -8.8979 10.3103 -0.4021
 -8.8503 10.2982 -0.4009
 -8.8934 10.3051 -0.3204
 -8.8922 10.3037 -0.3744
 -8.8957 10.3077 -0.3475
 -8.8903 10.3014 -0.2660
 -8.8403 10.2866 -0.3192
 -8.8439 10.2908 -0.3733
 -8.8447 10.2917 -0.3463
 -8.7987 10.2798 -0.3454
 -8.7972 10.2781 -0.3185
 -8.8420 10.2886 -0.2651
 -8.8414 10.2879 -0.2920
 -8.7972 10.2781 -0.2914
 -8.8914 10.3028 -0.2118
 -8.8913 10.3026 -0.2390
 -8.8911 10.3024 -0.1848
 -8.8431 10.2899 -0.2381
];
 
plane3 = [
 -8.7668 10.3218 -0.4267
 -8.8744 10.3263 -0.4560
 -8.7807 10.3382 -0.4543
 -8.8118 10.3373 -0.4551
 -8.8147 10.2984 -0.4003
 -8.8596 10.3091 -0.4282
 -8.8263 10.3121 -0.4275
 -8.7545 10.3073 -0.3992
 -8.7821 10.3024 -0.3996
 -8.7932 10.3154 -0.4269
 -8.7576 10.2737 -0.3445
 -8.8032 10.2850 -0.3726
 -8.7698 10.2880 -0.3719
 -8.6894 10.3576 -0.4529
 -8.7478 10.3424 -0.4536
 -8.7209 10.3506 -0.4534
 -8.7354 10.3278 -0.4262
 -8.7101 10.3377 -0.4259
 -8.7437 10.2946 -0.3717
];
 
plane4 = [
 -8.6503 10.3510 -0.3985
 -8.6315 10.3681 -0.3985
 -8.6334 10.3705 -0.4520
 -8.6954 10.3203 -0.3985
 -8.6694 10.3337 -0.3982
 -8.6744 10.3396 -0.4252
 -8.6550 10.3567 -0.4252
 -8.6332 10.3702 -0.4251
 -8.7059 10.2929 -0.3174
 -8.7298 10.2782 -0.3441
 -8.7117 10.2996 -0.3713
 -8.7069 10.2940 -0.3442
 -8.6881 10.3115 -0.3442
 -8.7283 10.2765 -0.3173
 -8.6858 10.3089 -0.3174
 -8.6907 10.3147 -0.3712
 -8.6435 10.3430 -0.3173
 -8.6684 10.3324 -0.3711
 -8.6662 10.3298 -0.3440
 -8.6463 10.3463 -0.3711
 -8.6288 10.3650 -0.3712
 -8.6454 10.3453 -0.3443
 -8.6273 10.3631 -0.3443
 -8.7063 10.2933 -0.2904
 -8.6863 10.3095 -0.2636
 -8.6875 10.3109 -0.2905
 -8.6632 10.3263 -0.2636
 -8.6636 10.3268 -0.2903
 -8.6423 10.3415 -0.2903
 -8.6206 10.3550 -0.2634
 -8.6220 10.3567 -0.3174
 -8.6235 10.3586 -0.2903
 -8.6664 10.3301 -0.2369
 -8.6401 10.3389 -0.2366
 -8.6232 10.3582 -0.2099
 -8.6199 10.3542 -0.2366    
];
 
% A-priori Coordinates of Apex
xapex = -8.7;
yapex = 10.3;
zapex = -0.3;
 
% Plane Parameters of face 1...4
nx1 = -0.211;
ny1 =  0.875;
nz1 = -0.436;
d1  = 10.971;
 
nx2 =  0.280;
ny2 =  0.959;
nz2 =  0.029;
d2  =  7.389;
 
nx3 = -0.152;
ny3 =  0.831;
nz3 =  0.536;
d3  =  9.678;
 
nx4 = -0.627;
ny4 =  0.774;
nz4 =  0.091;
d4  = 13.396;
 
% Vector of unknows
X0 = [xapex, yapex, zapex,...
      nx1, ny1, nz1, d1,...
      nx2, ny2, nz2, d2,...
      nx3, ny3, nz3, d3,...
      nx4, ny4, nz4, d4]';
 
Qll = eye(3*(length(plane1)+length(plane2)+length(plane3)+length(plane4)));
 
plane10 = plane1;
plane20 = plane2;
plane30 = plane3;
plane40 = plane4;
 
maxIteration = 200;
dx = inf;
v = zeros(3*(length(plane1)+length(plane2)+length(plane3)+length(plane4)),1);
while (abs(max(dx)) > sqrt(eps) && maxIteration >= 0)
    maxIteration = maxIteration - 1;
 
    % Jacobi
    A = zeros(length(plane1)+length(plane2)+length(plane3)+length(plane4), length(X0));
    w = zeros(size(A,1), 1);
    % Restriction
    R = zeros(8, length(X0));
    r = zeros(8, 1);
    % Condition
    B = zeros(size(A,1), size(A,1)*3);
 
 
    A(1:length(plane1), 4:7)                                                                                             = [plane1 -ones(length(plane1), 1)];
    A(length(plane1)+1:length(plane1)+length(plane2), 8:11)                                                              = [plane2 -ones(length(plane2), 1)];
    A(length(plane1)+length(plane2)+1:length(plane1)+length(plane2)+length(plane3), 12:15)                               = [plane3 -ones(length(plane3), 1)];
    A(length(plane1)+length(plane2)+length(plane3)+1:length(plane1)+length(plane2)+length(plane3)+length(plane4), 16:19) = [plane4 -ones(length(plane4), 1)];
 
    w(1:length(plane1), 1)                                                                                           = plane1 * X0(4:6)   - X0(7);
    w(length(plane1)+1:length(plane1)+length(plane2), 1)                                                             = plane2 * X0(8:10)  - X0(11);
    w(length(plane1)+length(plane2)+1:length(plane1)+length(plane2)+length(plane3), 1)                               = plane3 * X0(12:14) - X0(15);
    w(length(plane1)+length(plane2)+length(plane3)+1:length(plane1)+length(plane2)+length(plane3)+length(plane4), 1) = plane4 * X0(16:18) - X0(19);
 
    B(1:length(plane1), 1:3*length(plane1)) = kron(eye(length(plane1)), X0(4:6)');
    B(length(plane1)+1:length(plane1)+length(plane2), 3*length(plane1)+1:3*(length(plane1)+length(plane2))) = kron(eye(length(plane2)), X0(8:10)');
    B(length(plane1)+length(plane2)+1:length(plane1)+length(plane2)+length(plane3), 3*(length(plane1)+length(plane2))+1:3*(length(plane1)+length(plane2)+length(plane3))) = kron(eye(length(plane3)), X0(12:14)');
    B(length(plane1)+length(plane2)+length(plane3)+1:length(plane1)+length(plane2)+length(plane3)+length(plane4), 3*(length(plane1)+length(plane2)+length(plane3))+1:3*(length(plane1)+length(plane2)+length(plane3)+length(plane4))) = kron(eye(length(plane4)), X0(16:18)');
 
    R(1, 4:6)   = 2.*X0(4:6);
    R(2, 8:10)  = 2.*X0(8:10);
    R(3, 12:14) = 2.*X0(12:14);
    R(4, 16:18) = 2.*X0(16:18);
 
    R(5, [1:3,  4:7])  = [X0(4:6)'   X0(1:3)' -1];
    R(6, [1:3, 8:11])  = [X0(8:10)'  X0(1:3)' -1];
    R(7, [1:3,12:15])  = [X0(12:14)' X0(1:3)' -1];
    R(8, [1:3,16:19])  = [X0(16:18)' X0(1:3)' -1];
 
    r(1,1) = X0(4:6)'  *X0(4:6)   - 1;
    r(2,1) = X0(8:10)' *X0(8:10)  - 1;
    r(3,1) = X0(12:14)'*X0(12:14) - 1;
    r(4,1) = X0(16:18)'*X0(16:18) - 1;
 
    r(5,1) = X0(4:6)'  *X0(1:3) - X0(7);
    r(6,1) = X0(8:10)' *X0(1:3) - X0(11);
    r(7,1) = X0(12:14)'*X0(1:3) - X0(15);
    r(8,1) = X0(16:18)'*X0(1:3) - X0(19);
 
    w = -B*v + w;
 
    N = [   B*Qll*B'                          A           zeros(size(A,1),size(R,1))
               A'                      zeros(size(A,2))             R'
         zeros(size(R,1), size(A,1))      R           zeros(size(R,1),size(R,1))     ];
 
    invN = -inv(N);
    kxk = invN * [w; zeros(size(X0,1),1); r];
 
    k1 = kxk(1:1:size(A,1));
    dx = kxk(size(k1,1)+1:1:size(X0,1)+size(k1,1));
    k2 = kxk(size(k1,1)+size(dx,1)+1:1:size(kxk,1));
 
    X0 = X0+dx;
    v = Qll*B'*k1;
 
    for i=1:length(plane1)
        plane1(i, :) = plane10(i, :) + v(3*(i-1)+1:3*(i-1)+3)';
    end
 
    for i=1:length(plane2)
        plane2(i, :) = plane20(i, :) + v(3*length(plane1)+3*(i-1)+1 : 3*length(plane1)+3*(i-1)+3)';
    end
 
    for i=1:length(plane3)
        plane3(i, :) = plane30(i, :) + v(3*(length(plane1)+length(plane2))+3*(i-1)+1 : 3*(length(plane1)+length(plane2))+3*(i-1)+3)';
    end
 
    for i=1:length(plane4)
        plane4(i, :) = plane40(i, :) + v(3*(length(plane1)+length(plane2)+length(plane3))+3*(i-1)+1 : 3*(length(plane1)+length(plane2)+length(plane3))+3*(i-1)+3)';
    end
end
 
omega = v'*inv(Qll)*v + 2*k1'*(B*v+A*dx+w) + 2*k2'*(R*dx+r);
Qxx = invN(size(A,1)+1:1:size(A,1)+size(X0,1), size(A,1)+1:1:size(A,1)+size(X0,1));
sigma = omega / (size(A,1)+size(R,1)-size(A,2));
Qxx = sigma .* Qxx;  
 
disp([X0 sqrt(diag(Qxx))]);
 
plot3(plane1(:,1), plane1(:,2), plane1(:,3),'.b',...
      plane2(:,1), plane2(:,2), plane2(:,3),'.r',...
      plane3(:,1), plane3(:,2), plane3(:,3),'.g',...
      plane4(:,1), plane4(:,2), plane4(:,3),'.m',...
      X0(1), X0(2), X0(3), '*c');

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Tags:
Orthogonal Distance Fit, Formanalysis, Pyramid, Octave, Gauß-Helmert

Eintrag gesperrt

16516 Views

Hi, thanks a lot for your valuable effort it is really amazing algorithm. I run it is working very well. I will use this algorithm for my others pyramid that I collected from different point densities. So, could you please send to me your full name with the title of your algorithm, in order to use it as a reference in my thesis. for example, (Losle M., 2012), and the full reference needs to be written in full address in the page of references. What would you like to write to you in the reference page?
I will try to understand every element that you used in this algorithm.
Thanks again and happy Christmas and NEW YEAR.
Kindly regards

Eintrag gesperrt

15916 Views

Hi Balata,

So, could you please send to me your full name with the title of your algorithm, in order to use it as a reference in my thesis.

A suggestion is maybe:

Lösler, M. (2012). Orthogonal Distance Fit of a Pyramid. http://forum.diegeodaeten.de/index.php?id=4083 (last access: 2012-12-02).

I will try to understand every element that you used in this algorithm.

The used least squares algorithm is a common way for data fitting. It can also be used for other problems like transformations. It is often called Gauß-Helmert model or general least squares. An understandable description of the Gauß-Helmert model (in english) is given by e.g. Neitzel's Generalization of total least-squares on example of unweighted and weighted 2D similarity transformation (Springer allows free download of PDFs til the end of this year!). Moreover, the work of Lenzmann & Lenzmann is to point out: Strenge Auswertung des nichtlinearen Gauß-Helmert-Modells.

To fit the pyramid, the described algorithm fits four planes in one pour. These four planes are parameterized by the so-called Hesse normal form. The first 4 restrictions normalize the normal vector of each plane (subsequent shown for plane 1).

R(1, 4:6)   = 2.*X0(4:6); % nx*nx + ny*ny + nz*nz == 1
...
r(1,1) = X0(4:6)'  *X0(4:6)   - 1;

To estimate the apex within the data fit, 4 additional restrictions are needed. These restrictions simple define the apex as a point of each plane (subsequent shown for plane 1).

R(5, [1:3,  4:7])  = [X0(4:6)'   X0(1:3)' -1]; % nx*xApex + ny*yApex + nz*zApex == d
...
r(5,1) = X0(4:6)'  *X0(1:3) - X0(7);

Good luck with your thesis. Maybe, you can forward a digital copy of your work after finishing to me.

Kind Regards,
Micha

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Tags:
Orthogonal Distance Fit, Pyramid

Eintrag gesperrt

16283 Views

Hi Meche,

Thanks a lot for all your valuable information and guides. I'm now is about to finish my thesis I have to submit my thesis by January- February 2013. surely, once, I finish my thesis, I will send a copy of it to you.
The problem of the precision of the apex was remained until I contacted with you, and you gave me a good indication of how precisely estimate. I did a program through matlab code using plane equation to determine the apex. The result of apex coordinates is quite similar to your result. But I couldn't determine the precision of the apex. A copy of the second part of the program(computation of the apex) I attached below.

I also would like to let you know for other targets such as cone and sphere I used JAG3D, as another method for fitting cone and sphere. Is this software is also refer to your self? if yes, can you write the source as you did for pyramid. This software I used since February 2012.

Can I ask if possible to send to me the program of JAG3D for fitting cone and sphere instead of software in order to be able to use through matlab. As the same as you did for pyramid.

Many thanks indeed. You are really a good person I saw it yet during the period of my PhD study.

THEINTERSECTIONPOINT =
 
   -8.7482
   10.2641
   -0.3286
 
 clear all
clc
C1=load('D:\New Folder\pos114.10.11\p11x1py4\c1p1x1.txt');
X=C1(:,1);
Y=C1(:,2);
Z=C1(:,3);
figure(1)
plot3(X,Y,Z,('k.'))
xlabel('x-direction'),ylabel('y-direction'),zlabel('z-direction')
grid 
 
x = C1(:,1);
y = C1(:,2);
a=-0.0193;b=0.0797;c=-0.0396;
z = (1-a*x-b*y)/c; % Equation of the plane containing 
% (x,y,z) points is a*x+b*y+c*z=1
 
Xcolv = x(:); % Make X a column vector
Ycolv = y(:); % Make Y a column vector
Zcolv = z(:); % Make Z a column vector
Const = ones(size(Xcolv)); % Vector of ones for constant term
 
Coefficients = [Xcolv Ycolv Const]\Zcolv; % Find the coefficients
XCoeff = Coefficients(1); % X coefficient
YCoeff = Coefficients(2); % X coefficient
CCoeff = Coefficients(3); % constant term
% Using the above variables, z = XCoeff * x + YCoeff * y + CCoeff
figure(2)
L=plot3(x,y,z,'ko'); % Plot the original data points
set(L,'Markersize',0.5*get(L,'Markersize')) % Making the circle markers larger
set(L,'Markerfacecolor','k') % Filling in the markers
hold on
grid
[xx, yy]=meshgrid(min(X):0.01:max(X),min(Y):0.01:max(Y)); % Generating a regular grid for plotting
zz = XCoeff * xx + YCoeff * yy + CCoeff;
d1=z-C1(:,3);
sdvsurf1=std(zz);
sdv1=std(d1)
surf(xx,yy,zz) % Plotting the surface
title(sprintf('Plotting plane z=(%f)*x+(%f)*y+(%f)',XCoeff, YCoeff, CCoeff))
 
D1=load('D:\New Folder\pos114.10.11\p11x1py4\d1p1x1.txt');
X=D1(:,1);
Y=D1(:,2);
Z=D1(:,3);
figure(3)
plot3(X,Y,Z,('r.'))
xlabel('x-direction'),ylabel('y-direction'),zlabel('z-direction')
grid on
x = D1(:,1);
y = D1(:,2);
a=0.0373;b=0.1293;c=0.0040;
z = (1-a*x-b*y)/c; % Equation of the plane containing 
% (x,y,z) points is a*x+b*y+c*z=1
 
Xcolv = x(:); % Make X a column vector
Ycolv = y(:); % Make Y a column vector
Zcolv = z(:); % Make Z a column vector
Const = ones(size(Xcolv)); % Vector of ones for constant term
 
Coefficients = [Xcolv Ycolv Const]\Zcolv; % Find the coefficients
XCoeff = Coefficients(1); % X coefficient
YCoeff = Coefficients(2); % X coefficient
CCoeff = Coefficients(3); % constant term
% Using the above variables, z = XCoeff * x + YCoeff * y + CCoeff
figure(4)
L=plot3(x,y,z,'ro'); % Plot the original data points
set(L,'Markersize',0.5*get(L,'Markersize')) % Making the circle markers larger
%set(L,'Markerfacecolor','r') % Filling in the markers
hold on
grid
[xx, yy]=meshgrid(min(X):0.01:max(X),min(Y):0.01:max(Y)); % Generating a regular grid for plotting
zz = XCoeff * xx + YCoeff * yy + CCoeff;
d2=z-D1(:,3);
sdvsur2=std(zz);
sdv2=std(d2)
surf(xx,yy,zz) % Plotting the surface
title(sprintf('Plotting plane z=(%f)*x+(%f)*y+(%f)',XCoeff, YCoeff, CCoeff))
 
E1=load('D:\New Folder\pos114.10.11\p11x1py4\e1p1x1.txt');
X=E1(:,1);
Y=E1(:,2);
Z=E1(:,3);
figure(5)
plot3(X,Y,Z,('b.'))
xlabel('x-direction'),ylabel('y-direction'),zlabel('z-direction')
grid on
x = E1(:,1);
y = E1(:,2);
a=-0.0468;b=0.0576;c=0.0068;
z = (1-a*x-b*y)/c; % Equation of the plane containing 
% (x,y,z) points is a*x+b*y+c*z=1
 
Xcolv = x(:); % Make X a column vector
Ycolv = y(:); % Make Y a column vector
Zcolv = z(:); % Make Z a column vector
Const = ones(size(Xcolv)); % Vector of ones for constant term
 
Coefficients = [Xcolv Ycolv Const]\Zcolv; % Find the coefficients
XCoeff = Coefficients(1); % X coefficient
YCoeff = Coefficients(2); % X coefficient
CCoeff = Coefficients(3); % constant term
% Using the above variables, z = XCoeff * x + YCoeff * y + CCoeff
figure(6)
L=plot3(x,y,z,'b*'); % Plot the original data points
set(L,'Markersize',0.5*get(L,'Markersize')) % Making the circle markers larger
set(L,'Markerfacecolor','b') % Filling in the markers
hold on
grid
[xx, yy]=meshgrid (min(X):0.01:max(X),min(Y):0.01:max(Y),min(Z):0.01:max(Z)); % Generating a regular grid for plotting
zz = XCoeff * xx + YCoeff * yy + CCoeff;
d3=z-E1(:,3);
sdvsur3=std(zz);
sdv3=std(d3)
surf(xx,yy,zz) % Plotting the surface
title(sprintf('Plotting plane z=(%f)*x+(%f)*y+(%f)',XCoeff, YCoeff, CCoeff))
% By rotating the surface, you can see that the points lie on the plane
% Also, if we multiply both sides of the equation in the title by 4,
% we get the equation in the comment on the third line of this example
 
F1=load('D:\New Folder\pos114.10.11\p11x1py4\f1p1x1.txt');
X=F1(:,1);
Y=F1(:,2);
Z=F1(:,3);
figure(7)
plot3(X,Y,Z,('g.'))
xlabel('x-direction'),ylabel('y-direction'),zlabel('z-direction')
grid on
x = F1(:,1);
y = F1(:,2);
a=-0.0158;b=0.0857;c=0.0551;
z = (1-a*x-b*y)/c; % Equation of the plane containing 
% (x,y,z) points is a*x+b*y+c*z=1
 
Xcolv = x(:); % Make X a column vector
Ycolv = y(:); % Make Y a column vector
Zcolv = z(:); % Make Z a column vector
Const = ones(size(Xcolv)); % Vector of ones for constant term
 
Coefficients = [Xcolv Ycolv Const]\Zcolv; % Find the coefficients
XCoeff = Coefficients(1); % X coefficient
YCoeff = Coefficients(2); % X coefficient
CCoeff = Coefficients(3); % constant term
% Using the above variables, z = XCoeff * x + YCoeff * y + CCoeff
figure(8)
L=plot3(x,y,z,'g+'); % Plot the original data points
set(L,'Markersize',0.5*get(L,'Markersize')) % Making the circle markers larger
%set(L,'Markerfacecolor','g') % Filling in the markers
hold on
grid
[xx, yy]=meshgrid(min(X):0.01:max(X),min(Y):0.01:max(Y)); % Generating a regular grid for plotting
zz = XCoeff * xx + YCoeff * yy + CCoeff;
d4=z-F1(:,3);
sdvsur4=std(zz);
sdv4=std(d4)
surf(xx,yy,zz) % Plotting the surface
title(sprintf('Plotting plane z=(%f)*x+(%f)*y+(%f)',XCoeff, YCoeff, CCoeff))
sdvapex=sqrt(sdv1^2+sdv2^2+sdv3^2+sdv4)/4
 
C1=load('D:\New Folder\pos114.10.11\p11x1py4\c1p1x1.txt');
X=C1(:,1);
Y=C1(:,2);
Z=C1(:,3);
figure(9)
plot3(X,Y,Z,('k.'))
xlabel('x-direction'),ylabel('y-direction'),zlabel('z-direction')
grid 
hold on
 
D1=load('D:\New Folder\pos114.10.11\p11x1py4\d1p1x1.txt');
X=D1(:,1);
Y=D1(:,2);
Z=D1(:,3);
plot3(X,Y,Z,('ro'))
 
E1=load('D:\New Folder\pos114.10.11\p11x1py4\e1p1x1.txt');
X=E1(:,1);
Y=E1(:,2);
Z=E1(:,3);
plot3(X,Y,Z,('b*'))
 
F1=load('D:\New Folder\pos114.10.11\p11x1py4\f1p1x1.txt');
X=F1(:,1);
Y=F1(:,2);
Z=F1(:,3);
plot3(X,Y,Z,('g+'))
b=ones(length(C1),1);
par1=inv(C1'*C1)*C1'*b % first parameters
standreddivationofsurfaceC=std(C1);
%the dark grey coulour, normal tape
 
b=ones(length(D1),1);
par2=inv(D1'*D1)*D1'*b  %second parameters
standreddivationofsurfaceD=std(D1);
 
%the light grey colour, high reflectivity tape
b=ones(length(E1),1);
par3=inv(E1'*E1)*E1'*b  %third parameters
standreddivationofsurfaceE=std(E1);
%the green coulour, normal tape
 
b=ones(length(F1),1);
par4=inv(F1'*F1)*F1'*b  %second parameters
standreddivationofsurfaceF1=std(F1)
Q4=[par1';par2';par3';par4'];
 
%THE INTERSECTION POINT
J=[1;1;1;1];
THEINTERSECTIONPOINT =inv(Q4'*Q4)*Q4'*J
T=THEINTERSECTIONPOINT;
plot3(T(1,1),T(2,1),T(3,1),'y+')

Eintrag gesperrt

16263 Views

Hi Balata,

surely, once, I finish my thesis, I will send a copy of it to you.

Thank you - don't forget it! ;-)

and you gave me a good indication of how precisely estimate.

I was surprised, how good the algorithm estimates the pyramid-surface and an uncertainty of about 2 mm seems to be realistic.

But I couldn't determine the precision of the apex. A copy of the second part of the program(computation of the apex) I attached below.

I don't study the complete code but it seems, that you estimate each plane in a single step and combine these solutions to derive the apex. The "problem" is the imperfection of the measurement. Therefore, it is possible that no intersection point is given by the 3 (or more) planes. Thus, the consideration was to formulate the apex as restriction and forcing an intersection. So, the plane parameters maybe deviate from a single solution.

This procedure is not new or developed by me. A similar problem was solved by Hermann Bähr: Präzise Vermessung des Phasenreferenzpunktes von Corner-Reflektoren pp. 13-24 (link to PDF).

Is [JAG3D] is also refer to your self?

I'm the developer, yes. I don't have a publication which focus the software itself because it is a leisure project. I introduce JAG3D at the last tech12-continuing education seminar:

Lösler, M. (2012). JAG3D – Ein kostenfreies Programm zur Netzausgleichung und
Deformationsanalyse. Neitzel, F. (Hrsg.): tech12 - Aktuelle Trends der Messdatenauswertung
in Kataster- und Ingenieurvermessung, 19./20. April 2012, Berlin, Deutschland.

But you can also write in a similar way to the pyramid:

Lösler, M. (2012). JAG3D – A free Program for Network Adjustment and Deformation Analysis.
http://derletztekick.com/software/netzausgleichung (last access: 2012-12-02).

Can I ask if possible to send to me the program of JAG3D for fitting cone and sphere instead of software in order to be able to use through matlab.

The source code is free. You can find it e.g. at the SVN on sf.net.

Together with Dr. Matrin Nitschke, we publish some Matlab/Octave code to fit a plane or ellipse but I never create a m-File for cone or sphere.

kind regards
Micha

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Eintrag gesperrt

16182 Views

Hi Meche,

Many thanks for your effort and help.

I used your algorithm for my pyramid targets (4 faces). The algorithm has worked very well. However, when it applied to 3 faces after removing one plane an error occur in N matrix. As a result, I am still checking each step to see what the problem is with the 3 faces.
The algorithm is little bit hard to understand therefore I'm not an expert in this field,

Also, if I have trouble understanding your algorithm, would I be able to ask you a question?

Regards

Eintrag gesperrt

15968 Views

Hi,

Also, if I have trouble understanding your algorithm, would I be able to ask you a question?

Surely, I will answered at the weekend.

Only short, I'm in a hurry. I just remove the plane 4, the modified code is:

%% 3-Faces-Pyramid
%
%   nx(i) * x(j,i) + ny(i) * y(j,i) + nz(i) * z(j,i) - d(i) = 0
% with i = 1 .. 3
%
%   ||n(i)|| == 1
% 
%  and 
%
%   n(i)' * Apex -d(i) == 0 
%
%
% @version 1.0
% @author Michael Loesler - derletztekick.com
 
clear all
close all
format long g
clc
 
plane1 = [
 -8.7740 10.2929 -0.2645
 -8.7546 10.2701 -0.3177
 -8.8042 10.2864 -0.2647
 -8.7602 10.2767 -0.2908
 -8.7154 10.3041 -0.2639
 -8.7347 10.2841 -0.2905
 -8.7470 10.2985 -0.2641
 -8.8216 10.3067 -0.2113
 -8.8498 10.2977 -0.2114
 -8.7888 10.3102 -0.2109
 -8.8108 10.2940 -0.2379
 -8.7844 10.3051 -0.2378
 -8.7524 10.3050 -0.2374
 -8.7637 10.3184 -0.2109
 -8.8621 10.3120 -0.1847
 -8.8280 10.3142 -0.1843
 -8.7998 10.3231 -0.1842
 -8.7734 10.3298 -0.1841
 -8.7055 10.3322 -0.2103
 -8.7251 10.3155 -0.2372
 -8.6989 10.3245 -0.2371
 -8.7321 10.3238 -0.2106
 -8.6751 10.3405 -0.2103
 -8.7436 10.3375 -0.1838
 -8.7142 10.3427 -0.1837
 -8.6877 10.3555 -0.1837
 -8.6584 10.3609 -0.1835
];
 
plane2 = [
 -8.9011 10.3140 -0.4292
 -8.8979 10.3103 -0.4021
 -8.8503 10.2982 -0.4009
 -8.8934 10.3051 -0.3204
 -8.8922 10.3037 -0.3744
 -8.8957 10.3077 -0.3475
 -8.8903 10.3014 -0.2660
 -8.8403 10.2866 -0.3192
 -8.8439 10.2908 -0.3733
 -8.8447 10.2917 -0.3463
 -8.7987 10.2798 -0.3454
 -8.7972 10.2781 -0.3185
 -8.8420 10.2886 -0.2651
 -8.8414 10.2879 -0.2920
 -8.7972 10.2781 -0.2914
 -8.8914 10.3028 -0.2118
 -8.8913 10.3026 -0.2390
 -8.8911 10.3024 -0.1848
 -8.8431 10.2899 -0.2381
];
 
plane3 = [
 -8.7668 10.3218 -0.4267
 -8.8744 10.3263 -0.4560
 -8.7807 10.3382 -0.4543
 -8.8118 10.3373 -0.4551
 -8.8147 10.2984 -0.4003
 -8.8596 10.3091 -0.4282
 -8.8263 10.3121 -0.4275
 -8.7545 10.3073 -0.3992
 -8.7821 10.3024 -0.3996
 -8.7932 10.3154 -0.4269
 -8.7576 10.2737 -0.3445
 -8.8032 10.2850 -0.3726
 -8.7698 10.2880 -0.3719
 -8.6894 10.3576 -0.4529
 -8.7478 10.3424 -0.4536
 -8.7209 10.3506 -0.4534
 -8.7354 10.3278 -0.4262
 -8.7101 10.3377 -0.4259
 -8.7437 10.2946 -0.3717
];
 
 
% A-priori Coordinates of Apex
xapex = -8.7;
yapex = 10.3;
zapex = -0.3;
 
% Plane Parameters of face 1...4
nx1 = -0.211;
ny1 =  0.875;
nz1 = -0.436;
d1  = 10.971;
 
nx2 =  0.280;
ny2 =  0.959;
nz2 =  0.029;
d2  =  7.389;
 
nx3 = -0.152;
ny3 =  0.831;
nz3 =  0.536;
d3  =  9.678;
 
 
% Vector of unknows
X0 = [xapex, yapex, zapex,...
      nx1, ny1, nz1, d1,...
      nx2, ny2, nz2, d2,...
      nx3, ny3, nz3, d3]';
 
Qll = eye(3*(length(plane1)+length(plane2)+length(plane3)));
 
plane10 = plane1;
plane20 = plane2;
plane30 = plane3;
 
maxIteration = 200;
dx = inf;
v = zeros(3*(length(plane1)+length(plane2)+length(plane3)),1);
while (abs(max(dx)) > sqrt(eps) && maxIteration >= 0)
    maxIteration = maxIteration - 1;
 
    % Jacobi
    A = zeros(length(plane1)+length(plane2)+length(plane3), length(X0));
    w = zeros(size(A,1), 1);
    % Restriction
    R = zeros(6, length(X0));
    r = zeros(6, 1);
    % Condition
    B = zeros(size(A,1), size(A,1)*3);
 
 
    A(1:length(plane1), 4:7)                                                                                             = [plane1 -ones(length(plane1), 1)];
    A(length(plane1)+1:length(plane1)+length(plane2), 8:11)                                                              = [plane2 -ones(length(plane2), 1)];
    A(length(plane1)+length(plane2)+1:length(plane1)+length(plane2)+length(plane3), 12:15)                               = [plane3 -ones(length(plane3), 1)];
 
    w(1:length(plane1), 1)                                                                                           = plane1 * X0(4:6)   - X0(7);
    w(length(plane1)+1:length(plane1)+length(plane2), 1)                                                             = plane2 * X0(8:10)  - X0(11);
    w(length(plane1)+length(plane2)+1:length(plane1)+length(plane2)+length(plane3), 1)                               = plane3 * X0(12:14) - X0(15);
 
    B(1:length(plane1), 1:3*length(plane1)) = kron(eye(length(plane1)), X0(4:6)');
    B(length(plane1)+1:length(plane1)+length(plane2), 3*length(plane1)+1:3*(length(plane1)+length(plane2))) = kron(eye(length(plane2)), X0(8:10)');
    B(length(plane1)+length(plane2)+1:length(plane1)+length(plane2)+length(plane3), 3*(length(plane1)+length(plane2))+1:3*(length(plane1)+length(plane2)+length(plane3))) = kron(eye(length(plane3)), X0(12:14)');
 
    R(1, 4:6)   = 2.*X0(4:6);
    R(2, 8:10)  = 2.*X0(8:10);
    R(3, 12:14) = 2.*X0(12:14);
 
    R(4, [1:3,  4:7])  = [X0(4:6)'   X0(1:3)' -1];
    R(5, [1:3, 8:11])  = [X0(8:10)'  X0(1:3)' -1];
    R(6, [1:3,12:15])  = [X0(12:14)' X0(1:3)' -1];
 
    r(1,1) = X0(4:6)'  *X0(4:6)   - 1;
    r(2,1) = X0(8:10)' *X0(8:10)  - 1;
    r(3,1) = X0(12:14)'*X0(12:14) - 1;
 
    r(4,1) = X0(4:6)'  *X0(1:3) - X0(7);
    r(5,1) = X0(8:10)' *X0(1:3) - X0(11);
    r(6,1) = X0(12:14)'*X0(1:3) - X0(15);
 
    w = -B*v + w;
 
    N = [   B*Qll*B'                          A           zeros(size(A,1),size(R,1))
               A'                      zeros(size(A,2))             R'
         zeros(size(R,1), size(A,1))      R           zeros(size(R,1),size(R,1))     ];
 
    invN = -inv(N);
    kxk = invN * [w; zeros(size(X0,1),1); r];
 
    k1 = kxk(1:1:size(A,1));
    dx = kxk(size(k1,1)+1:1:size(X0,1)+size(k1,1));
    k2 = kxk(size(k1,1)+size(dx,1)+1:1:size(kxk,1));
 
    X0 = X0+dx;
    v = Qll*B'*k1;
 
    for i=1:length(plane1)
        plane1(i, :) = plane10(i, :) + v(3*(i-1)+1:3*(i-1)+3)';
    end
 
    for i=1:length(plane2)
        plane2(i, :) = plane20(i, :) + v(3*length(plane1)+3*(i-1)+1 : 3*length(plane1)+3*(i-1)+3)';
    end
 
    for i=1:length(plane3)
        plane3(i, :) = plane30(i, :) + v(3*(length(plane1)+length(plane2))+3*(i-1)+1 : 3*(length(plane1)+length(plane2))+3*(i-1)+3)';
    end
end
 
omega = v'*inv(Qll)*v + 2*k1'*(B*v+A*dx+w) + 2*k2'*(R*dx+r);
Qxx = invN(size(A,1)+1:1:size(A,1)+size(X0,1), size(A,1)+1:1:size(A,1)+size(X0,1));
sigma = omega / (size(A,1)+size(R,1)-size(A,2));
Qxx = sigma .* Qxx;  
 
disp([X0 sqrt(diag(Qxx))]);
 
plot3(plane1(:,1), plane1(:,2), plane1(:,3),'.b',...
      plane2(:,1), plane2(:,2), plane2(:,3),'.r',...
      plane3(:,1), plane3(:,2), plane3(:,3),'.g',...
      X0(1), X0(2), X0(3), '*c');

regards
Micha

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Eintrag gesperrt

17476 Views

Hi Micha,

Thanks a lot, the program is working very well.
Can I ask what do you mean by "restriction"? For each plane you used two restrictions. I have gone through the Hesse normal form to understand. But, unfortunately I didn't get the point.
Regards

Eintrag gesperrt

15934 Views

Hi,

Can I ask what do you mean by "restriction"? For each plane you used two restrictions.

The first restriction normalizes the normal vector n of the plane. If you don't normalize the vector, the estimation will failed, because n is not unique. A vector [1 2 3]' is equal to [2 4 6]'. The lengths are changed but the directions are equal to each other.

The second restriction force the apex into the plane. This restriction is specific to your pyramid-problem.

regards
Micha

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Eintrag gesperrt

15980 Views

Hi Micha,

sorry for being late on replying to your explanation about restriction, due to I was far away from my office for a while. I highly appreciated for responding all my questions.

Can I ask, if you could summarize or describe the algorithm "pyramid" in a logical steps, in order to understand it. Because, some steps still I couldn't understand it.
Many thanks

Eintrag gesperrt

15959 Views

Hi,

Can I ask, if you could summarize or describe the algorithm "pyramid" in a logical steps, in order to understand it. Because, some steps still I couldn't understand it.

You are referred to the listed publications for details to the least-squares algorithm.

To fit the pyramid, the algorithm estimates the parameters of the 4 planes. The planes are parametrized as Hesse normal form:

$n_xx_i + n_yy_i + n_zz_i = d$

with the secondary condition
$||n|| \stackrel{!}{=} \sqrt{n_x^2 + n_y^2 + n_z^2} \stackrel{!}{=} 1$

For each plane the unknowns are $n_x$ , $n_y$ , $n_z$ and $d$ . In your case, the algorithm estimates 4×4 = 16 plane-parameters. Up to this point, the algorithm produces the same results as a single fit of each plane.

To estimate the intersection, the number of unknowns increased by 3 - the coordinates of the intersection-point. To force a intersection of the 4 planes, 4 additional conditions are introduced. If $K$ is the intersection-point, each plane has to contain $K$ , or similar:

$n_xx_K + n_yy_K + n_zz_K \stackrel{!}{=} d$

This is the whole mathematic model, which can be used in a common orthogonal distance fit algorithm - see listed publications to get formulas of Gauß-Helmert model.

kind regads
Micha

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Eintrag gesperrt

16196 Views

Hi Micha,

Thanks for your nice response, I have got the point.
I confused with the word of omega in the last step of the algorithm. What is "omega" represents?.
Regards

Eintrag gesperrt

16001 Views

Hi Balata,

What is "omega" represents?

Ω represents the objective function - see Eq. 7.

regards
Micha

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Eintrag gesperrt

16022 Views

Hi Michael,
I haven't being contact you for a long time because I was so busy with my writing phD thesis. Finally, I have submitted my thesis. And as I told you after I finish the viva, I will send a final copy of it to you. I'm now preparing myself for the exam

I have a question related to the LM algorithm. As I used JAG3D v.3.3 for fitting my targets. As I knew there was a list of Levenberg- Marquardt damping value μ. the algorithm acts to the Guass Newton (GN) method for very small value of damping parameter, and acts as the method of gradient decent for a large value of damping. The question here what is the damping value? could you please help to explain this LM damping value, because LM algorithm is adaptive by controlling its damping value. Therefore, this factor is important to be asked.
Kindly regards

Eintrag gesperrt

15837 Views

Hi,

after I finish the viva, I will send a final copy of it to you.

Yes, I still waiting. ;-)

Maybe a good explanation of the LM is given in Numerical Recipes in C p.684 or at Winipedia.

The question here what is the damping value?

It is just a scaling parameter to switch over to the method of gradient decent for a large value or, if this value is small, to the Gauß-Newton method. As seen at page 684 in NR in C, the damping parameter is a user-defined value. If the first approximation is far away from the minimum, you will start with a large value. If you can start with an appropriate approximation, you don't need a (large) damping parameter (choose a very small one) and switch on the powerful GN-method.

regards
Micha

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Eintrag gesperrt

15691 Views

Hi Mechae,
Thanks a lot for your nice answer. I also have a problem with the software JAG3D v.3.3, it doesn't work, when run it the message will appear " could not start the application" even with the new version the same message. I don't know, the problem in my computer " operating system" or the in the software, which might be needs some files to be downloaded separately such as crack files.

In orthogonal distance fit pyramid you said " I used a combined algorithm which estimates the 4 plane parameters and the apex in one pour, formulated as a Gauß-Helmert model with unknowns in the restrictions. The Octave(quick and dirty) code"
what are the two algorithms? what is the Octave(quick and dirty) code (i.e. why its name quick and dirty).

Regards

Eintrag gesperrt

15828 Views

Hi,

it doesn't work, when run it the message will appear " could not start the application"

Do you have installed Java? Your need at least version 1.6 or higher.

If I open a shell to check my installation, I get:

C:\>java -version
java version "1.7.0_11"
Java(TM) SE Runtime Environment (build 1.7.0_11-b21)
Java HotSpot(TM) 64-Bit Server VM (build 23.6-b04, mixed mode)
 
C:\>

what are the two algorithms?

Combined algorithm does not implicate the word two. It means, that the algorithm fits four planes together with the apex of your pyramid. So, you don't have to estimate each single plane and the apex in a separated way - you get it from one source.

what is the Octave(quick and dirty) code (i.e. why its name quick and dirty).

Octave is a powerful application for numerical computations. The phrase quick and dirty refers to my uncomment code with hard-coded values etc. This code was just a test case in order to show you a possible strategy.

Have a nice weekend
Micha

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Tags:
JAG3D, Java, Pyramid, Code, Octave

Eintrag gesperrt

15998 Views

Hi Michae,

Thanks a lot, I have downloded the java, and now the JAG3D software is working. This was my mistake when I had formatted my computer, I didn't know there was a link between them.

Many thanks.

Eintrag gesperrt

15720 Views

Hi,

I didn't know there was a link between them.

No problem. By the way; the J in JAG3D means JAVA. ;-)

reagrds
Micha

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Eintrag gesperrt

15772 Views

Hi Michael,

I have wanted to let you know that I did my PhD viva voce exam on 17th June 2013, and I have successfully passed with minor corrections within three months. As I told you an electonic copy of the final thesis I will send it to you.
Many thanks

Eintrag gesperrt

15514 Views

Hi Dr. Balata ;-)

Congratulations on passing the final! I am pleased to be able to read your work.

kind regards!
Micha

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Eintrag gesperrt

15583 Views

Hi,

to avoid SPAM-problems: If you send your thesis, please write a short notice here in the forum.

thank you
Micha

--
applied-geodesy.org - OpenSource Least-Squares Adjustment Software for Geodetic Sciences

Eintrag gesperrt

15609 Views

Hi Miche,
Thanks a lots. I am about to finish my correction next week, and I will send an electronic pdf copyit to your email via dropbox. you will get a link to open it.
Regards

Eintrag gesperrt

15535 Views

Pointing to my previous email, I've just wanted to let you know the algorithm that I used for the pyramid using Matlab code for computing the apex. The parameters of (a,b,c) for each plane I calculated separately by another algorithm. But the way for calculating stdv. for apex was wrong. YOU solved my problem by your nice algorithm

clear all
clc
C1=load('D:\New Folder\pos114.10.11\p11x1py4\c1p1x1.txt'); % point of plane 1
X=C1(:,1);
Y=C1(:,2);
Z=C1(:,3);
figure(1)
plot3(X,Y,Z,('k.'))
xlabel('x-direction'),ylabel('y-direction'),zlabel('z-direction')
grid 
 
x = C1(:,1);
y = C1(:,2);
a=-0.21131;b=0.87480;c=-0.43598;
z = (1-a*x-b*y)/c; % Equation of the plane containing 
% (x,y,z) points is a*x+b*y+c*z=1
 
Xcolv = x(:); % Make X a column vector
Ycolv = y(:); % Make Y a column vector
Zcolv = z(:); % Make Z a column vector
Const = ones(size(Xcolv)); % Vector of ones for constant term
 
Coefficients = [Xcolv Ycolv Const]\Zcolv; % Find the coefficients
XCoeff = Coefficients(1); % X coefficient
YCoeff = Coefficients(2); % X coefficient
CCoeff = Coefficients(3); % constant term
% Using the above variables, z = XCoeff * x + YCoeff * y + CCoeff
figure(2)
L=plot3(x,y,z,'ko'); % Plot the original data points
set(L,'Markersize',0.5*get(L,'Markersize')) % Making the circle markers larger
set(L,'Markerfacecolor','k') % Filling in the markers
hold on
grid
[xx, yy]=meshgrid(min(X):0.01:max(X),min(Y):0.01:max(Y)); % Generating a regular grid for plotting
zz = XCoeff * xx + YCoeff * yy + CCoeff;
d1=z-C1(:,3);
sdvsurf1=std(zz);
sdv1=std(d1)
surf(xx,yy,zz) % Plotting the surface
title(sprintf('Plotting plane z=(%f)*x+(%f)*y+(%f)',XCoeff, YCoeff, CCoeff))
 
D1=load('D:\New Folder\pos114.10.11\p11x1py4\d1p1x1.txt'); % points of plane 2
X=D1(:,1);
Y=D1(:,2);
Z=D1(:,3);
figure(3)
plot3(X,Y,Z,('r.'))
xlabel('x-direction'),ylabel('y-direction'),zlabel('z-direction')
grid on
x = D1(:,1);
y = D1(:,2);
a=0.280048226488374;b=0.959535677;c=0.02939857;
z = (1-a*x-b*y)/c; % Equation of the plane containing 
% (x,y,z) points is a*x+b*y+c*z=1
 
Xcolv = x(:); % Make X a column vector
Ycolv = y(:); % Make Y a column vector
Zcolv = z(:); % Make Z a column vector
Const = ones(size(Xcolv)); % Vector of ones for constant term
 
Coefficients = [Xcolv Ycolv Const]\Zcolv; % Find the coefficients
XCoeff = Coefficients(1); % X coefficient
YCoeff = Coefficients(2); % X coefficient
CCoeff = Coefficients(3); % constant term
% Using the above variables, z = XCoeff * x + YCoeff * y + CCoeff
figure(4)
L=plot3(x,y,z,'ro'); % Plot the original data points
set(L,'Markersize',0.5*get(L,'Markersize')) % Making the circle markers larger
%set(L,'Markerfacecolor','r') % Filling in the markers
hold on
grid
[xx, yy]=meshgrid(min(X):0.01:max(X),min(Y):0.01:max(Y)); % Generating a regular grid for plotting
zz = XCoeff * xx + YCoeff * yy + CCoeff;
d2=z-D1(:,3);
sdvsur2=std(zz);
sdv2=std(d2)
surf(xx,yy,zz) % Plotting the surface
title(sprintf('Plotting plane z=(%f)*x+(%f)*y+(%f)',XCoeff, YCoeff, CCoeff))
 
F1=load('D:\New Folder\pos114.10.11\p11x1py4\f1p1x1.txt'); % points of plane 3
X=F1(:,1);
Y=F1(:,2);
Z=F1(:,3);
figure(7)
plot3(X,Y,Z,('g.'))
xlabel('x-direction'),ylabel('y-direction'),zlabel('z-direction')
grid on
x = F1(:,1);
y = F1(:,2);
a=-0.15202644;b=0.830514033;c=0.535849234;
z = (1-a*x-b*y)/c; % Equation of the plane containing 
% (x,y,z) points is a*x+b*y+c*z=1
 
Xcolv = x(:); % Make X a column vector
Ycolv = y(:); % Make Y a column vector
Zcolv = z(:); % Make Z a column vector
Const = ones(size(Xcolv)); % Vector of ones for constant term
 
Coefficients = [Xcolv Ycolv Const]\Zcolv; % Find the coefficients
XCoeff = Coefficients(1); % X coefficient
YCoeff = Coefficients(2); % X coefficient
CCoeff = Coefficients(3); % constant term
% Using the above variables, z = XCoeff * x + YCoeff * y + CCoeff
figure(8)
L=plot3(x,y,z,'g+'); % Plot the original data points
set(L,'Markersize',0.5*get(L,'Markersize')) % Making the circle markers larger
%set(L,'Markerfacecolor','g') % Filling in the markers
hold on
grid
[xx, yy]=meshgrid(min(X):0.01:max(X),min(Y):0.01:max(Y)); % Generating a regular grid for plotting
zz = XCoeff * xx + YCoeff * yy + CCoeff;
d4=z-F1(:,3);
sdvsur4=std(zz);
sdv4=std(d4)
surf(xx,yy,zz) % Plotting the surface
title(sprintf('Plotting plane z=(%f)*x+(%f)*y+(%f)',XCoeff, YCoeff, CCoeff))
 
E1=load('D:\New Folder\pos114.10.11\p11x1py4\e1p1x1.txt'); %point of plane 4
X=E1(:,1);
Y=E1(:,2);
Z=E1(:,3);
figure(5)
plot3(X,Y,Z,('b.'))
xlabel('x-direction'),ylabel('y-direction'),zlabel('z-direction')
grid on
x = E1(:,1);
y = E1(:,2);
a=-0.626881329;b=0.773774228;c=0.09106725;
z = (1-a*x-b*y)/c; % Equation of the plane containing 
% (x,y,z) points is a*x+b*y+c*z=1
 
Xcolv = x(:); % Make X a column vector
Ycolv = y(:); % Make Y a column vector
Zcolv = z(:); % Make Z a column vector
Const = ones(size(Xcolv)); % Vector of ones for constant term
 
Coefficients = [Xcolv Ycolv Const]\Zcolv; % Find the coefficients
XCoeff = Coefficients(1); % X coefficient
YCoeff = Coefficients(2); % X coefficient
CCoeff = Coefficients(3); % constant term
% Using the above variables, z = XCoeff * x + YCoeff * y + CCoeff
figure(6)
L=plot3(x,y,z,'b*'); % Plot the original data points
set(L,'Markersize',0.5*get(L,'Markersize')) % Making the circle markers larger
set(L,'Markerfacecolor','b') % Filling in the markers
hold on
grid
[xx, yy]=meshgrid (min(X):0.01:max(X),min(Y):0.01:max(Y),min(Z):0.01:max(Z)); % Generating a regular grid for plotting
zz = XCoeff * xx + YCoeff * yy + CCoeff;
d3=z-E1(:,3);
sdvsur3=std(zz);
sdv3=std(d3)
surf(xx,yy,zz) % Plotting the surface
title(sprintf('Plotting plane z=(%f)*x+(%f)*y+(%f)',XCoeff, YCoeff, CCoeff))
% By rotating the surface, you can see that the points lie on the plane
% Also, if we multiply both sides of the equation in the title by 4,
% we get the equation in the comment on the third line of this example
 
sdvapex=sqrt(sdv1^2+sdv2^2+sdv3^2+sdv4)/4
 
C1=load('D:\New Folder\pos114.10.11\p11x1py4\c1p1x1.txt');
X=C1(:,1);
Y=C1(:,2);
Z=C1(:,3);
figure(9)
plot3(X,Y,Z,('k.'))
xlabel('x-direction'),ylabel('y-direction'),zlabel('z-direction')
grid 
hold on
 
D1=load('D:\New Folder\pos114.10.11\p11x1py4\d1p1x1.txt');
X=D1(:,1);
Y=D1(:,2);
Z=D1(:,3);
plot3(X,Y,Z,('ro'))
 
E1=load('D:\New Folder\pos114.10.11\p11x1py4\e1p1x1.txt');
X=E1(:,1);
Y=E1(:,2);
Z=E1(:,3);
plot3(X,Y,Z,('b*'))
 
F1=load('D:\New Folder\pos114.10.11\p11x1py4\f1p1x1.txt');
X=F1(:,1);
Y=F1(:,2);
Z=F1(:,3);
plot3(X,Y,Z,('g+'))

Eintrag gesperrt

16396 Views

RSS-Feed dieser Diskussion

Das GEO-Forum
http://forum.diegeodaeten.de

Abbreviations in JAG3D's FormFittingToolbox (Geodäsie/Vermessung)