* Bug 16365 fixed: median(m,'r'|'c') was wrong after 5dc990
[scilab.git] / scilab / modules / optimization / demos / datafit / ellipse.dem.sce
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2018 - Samuel GOUGEON
3 //
4 // This file is released under the 3-clause BSD license. See COPYING-BSD.
5
6 function demo_datafit()
7     function g = Gap(p, Data)
8         // p: a, b, center xc, yc, tilt °
9         C = cosd(p(5));
10         S = -sind(p(5));
11         x = Data(1,:) - p(3);
12         y = Data(2,:) - p(4);
13         g = ((x*C + y*S )/p(1)).^2 + ((-x*S + y*C)/p(2)).^2 - 1;
14     endfunction
15
16     // Generating the data
17     // -------------------
18     // Actual parameters :
19     //    a   b   xc  yc  tilt
20     a = grand(1,1,"unf",1,10);
21     b = grand(1,1,"unf",a/6, a);
22     C = grand(1,2,"unf",-10,10);
23     tilt = grand(1,1,"unf", -85,85);
24     Pa = [a b C(1) C(2) tilt];
25     Np = 30;            // Number of data points
26     Theta_min = grand(1,1,"unf",-180,-90)
27     Theta = grand(1,Np, "unf",Theta_min, Theta_min+270);
28     // untilted centered noisy arc
29     x = Pa(1)*(cosd(Theta) + grand(Theta, "unf",-0.07, 0.07));
30     y = Pa(2)*(sind(Theta) + grand(Theta, "unf",-0.07, 0.07));
31     // Tilting and off-centering the arc
32     A = Pa(5);
33     C = cosd(A); S = sind(A);
34     xe =  C*x + y*S + Pa(3);
35     ye = -S*x + y*C + Pa(4);
36     Data = [xe ; ye];
37
38     // Retrieving parameters from noisy data
39     // -------------------------------------
40     // Initial estimation
41     ab0 = (strange(xe)+strange(ye))/2;
42     xc0 = mean(xe);
43     yc0 = mean(ye);
44     A0 = -atand(reglin(xe,ye));
45     P0 = [ab0 ab0/2 xc0 yc0 A0];
46     // Parameters bounds
47     Pinf = [ 1     0.2     xc0-ab0/2, yc0-ab0/2 -89];
48     Psup = [1.2*ab0 1.2*ab0  xc0+ab0/2, yc0+ab0/2  89];// Fitting
49
50     // FITTING
51     [P, dmin] = datafit(Gap, Data, 'b',Pinf,Psup, P0);
52
53     // P(1) is the longest axis: postprocessing the orientation angle
54     if P(1)<P(2) then
55         P([1 2]) = P([2 1]);
56         if P(5)>0
57             P(5) = 90 - P(5);
58         else
59             P(5) = 90 + P(5);
60         end
61     end
62     if (P(5)*Pa(5))<0 & abs(-P(5)/Pa(5)-1)<0.3 & P(5)>5 then
63         P(5) = -P(5)
64     end
65
66     // COMPUTING THE FITTING CURVE
67     Theta = 0:2:360;
68     x = P(1) * cosd(Theta);
69     y = P(2) * sind(Theta);
70     A = P(5);
71     [x, y] = (x*cosd(A)+y*sind(A)+P(3), -x*sind(A)+y*cosd(A)+P(4));
72
73     // COMPUTING THE TRUE ELLIPSE CURVE
74     xt = Pa(1) * cosd(Theta);
75     yt = Pa(2) * sind(Theta);
76     At = Pa(5);
77     [xt, yt] = (xt*cosd(At)+yt*sind(At)+Pa(3), -xt*sind(At)+yt*cosd(At)+Pa(4));
78
79     // DISPLAY
80     // -------
81     my_handle = scf(100001);
82     drawlater
83     clf(my_handle, "reset");
84     gcf().figure_name = _("datafit() random demo : ellipse fitting (unweighted)");
85     demo_viewCode(get_absolute_file_path("ellipse.dem.sce")+"ellipse.dem.sce");
86
87     plot(xt, yt, "g")   // True ellipse
88     plot(xe, ye, "+")   // Data
89     plot(x, y, "r")     // Fitting ellipse
90     isoview
91     legend([_("True elliptical law")
92             _("Data = true + noise dx|y = +/-7%")
93             _("Best ellipse fitting Data")], "in_lower_right");
94     // Parameters values:
95     Patxt = msprintf("%5.2f   %5.2f   %5.2f   %5.2f  %5.1f", Pa);
96     Ptxt  = msprintf("%5.2f   %5.2f   %5.2f   %5.2f  %5.1f", P);
97     xstring(P(3), P(4), ["" "    a           b        xc       yc     A°"
98     "True:" Patxt
99     "Fit:" Ptxt
100     "" "    dmin = "+ msprintf("%5.3f", dmin)]);
101     set(gce(),"text_box",[0 0], "text_box_mode", "centered");
102     title(_("datafit: Fitting on 30 random unweighted data points"),"fontsize",3)
103     drawnow
104 endfunction
105
106 demo_datafit();
107 clear demo_datafit;