dea22f097bacc97c70fb751b0e59e24320173d5b
[scilab.git] / scilab / modules / elementary_functions / macros / unwrap.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) - 2013 - Samuel GOUGEON
3 //
4 // Copyright (C) 2012 - 2016 - Scilab Enterprises
5 //
6 // This file is hereby licensed under the terms of the GNU GPL v2.0,
7 // pursuant to article 5.3.4 of the CeCILL v.2.1.
8 // This file was originally licensed under the terms of the CeCILL v2.1,
9 // and continues to be available under such terms.
10 // For more information, see the COPYING file which you should have received
11 // along with this program.
12
13 function [retval, K] = unwrap(a, varargin)
14     // [unwr, cuspPoints] = unwrap(a)        // row or column vector or matrix
15     // [unwr, cuspPoints] = unwrap(a, jump)  // default jump = 2*%pi
16     // jump=0 => = average local slope
17     // [unwr, cuspPoints] = unwrap(a, jump, direc)
18     //          direc="c": unwrap/unfold along columns
19     //          direc="r": along rows
20     //          else: along both directions (default)
21     // [unwr, cuspPoints] = unwrap(a, "unfold")
22
23     // Initializations
24     // --------------
25     [lhs,rhs] = argn(0);
26     transposed = %f
27     unfold = %f     // unfold (1D) instead of unwrap (1D|2D)
28     jump = 2*%pi
29     direc  = ""
30     retval = []
31     K = []
32
33     // EXAMPLES
34     // --------
35     if rhs==0 then
36         %_unwrap()      // display 1D examples
37         halt(_("Press return to display 2D examples"))
38         %_unwrap("2D")  // display 2D examples
39         return
40     end
41
42     // =================== CHECKING PARAMETERS ====================
43     if type(a)~=1 then
44         msg = _("%s: Wrong type for input argument #%d: Real expected.\n")
45         error(msprintf(msg, "unwrap",1))
46     end
47     if rhs>1
48         if typeof(varargin(1))=="string" then
49             if varargin(1)=="unfold" then
50                 unfold = %t
51             else
52                 msg = _("%s: Wrong value for input argument #%d: Must be in the set  {%s}.\n")
53                 error(msprintf(msg, "unwrap",2,"""unfold"""))
54             end
55         else
56             jump = varargin(1)
57             if type(jump)~=1 then
58                 msg = _("%s: Wrong type for input argument #%d: Real expected.\n")
59                 error(msprintf(msg, "unwrap",2))
60             else
61                 jump = abs(jump(1))
62             end
63         end
64         if ~unfold & rhs>=3 then
65             direc = varargin(2)
66             if typeof(direc)~="string" then
67                 msg = _("%s: Wrong type for input argument #%d: String expected.\n")
68                 error(msprintf(msg, "unwrap", 3))
69             else
70                 direc = direc(1)
71                 if and(direc~=["r" "c"]) then
72                     msg = _("%s: Wrong value for input argument #%d: Must be in the set  {%s}.\n")
73                     error(msprintf(msg, "unwrap", 3, """r"",""c"""))
74                 end
75             end
76         end
77     end
78
79     if direc=="c" | size(a,2)==1
80         a = a.'
81         transposed = %t
82     end
83
84     [m,n] = size(a)
85     if (n < 4)
86         if transposed then
87             msg = _("%s: Wrong size for input argument #%d: Must have at least 4 rows.\n")
88         else
89             msg = _("%s: Wrong size for input argument #%d: Must have at least 4 columns.\n")
90         end
91         error(msprintf(msg, "unwrap",1))
92     end
93     if unfold & m>1 then
94         msg = _("%s: Wrong size for input argument #%d: Vector expected.\n")
95         error(msprintf(msg, "unwrap",1))
96     end
97
98     // ============================ PROCESSING ==============================
99
100     // Columns #1 and #$ are duplicated:
101     rae = [a(:,1)  a  a(:,$)];       // => n+2 columns
102     // Derivative along rows (assuming equally spaced x values)
103     d = rae(:,2:$) - rae(:, 1:$-1);     // n+1 columns
104     clear rae
105
106     // ---------------------------  U N W R A P  ---------------------------
107     if ~unfold then
108         // jump (may be)  (n-1 columns)
109         ju = d(:, 2:$-1)
110         // average Local derivative (before/after jump) (n-1 columns) :
111         avL = (d(:, 1:$-2)+d(:, 3:$))/2
112         // where
113         wh = abs(ju)>4*abs(avL) & d(:, 1:$-2).*d(:, 3:$)>0
114         wh = abs(ju)>5*abs(avL)
115         K = find(wh)
116         d = d(:, 2:$-1)                  // n-1 columns
117         if jump~=0 then
118             d(K) = ju(K) - sign(ju(K)).*jump
119             // Cleaning wrongly inserted jumps :
120             k = find(abs(d(K)-mean(d))>5*stdev(d))
121             if ~isempty(k)
122                 d(K(k)) = d(K(k)) - sign(d(K(k))).*jump
123             end
124         else
125             d(K) = avL(K)
126         end
127
128         if and(size(a)~=1) & direc==""
129             // 2D case: levels all rows according to the unwrapped 1st column
130             p = unwrap(a(:,1),jump)
131         else
132             p  = a(:,1)
133         end
134         r = [p d]
135
136         // ---------------------------  U N F O L D  ---------------------------
137     else
138         // Local radius of curvature: calculation
139         Rc = ((1+d(:,2:$).^2).^1.5) ./ abs(d(:,2:$)-d(:,1:$-1))
140         K = find(Rc < mean(Rc)/30) // criterium to detect cusp points
141         [nR,nC] = size(Rc)
142
143         // trimming edges (due to duplication of 1st & last columns):
144         tmp = find(K==1 | K==nC)
145         K(tmp) = []
146
147         // trimming duplicates (if any)
148         tmp = find(K(2:$)-K(1:$-1)==1)
149         K(tmp+1) = []
150
151         // fine targetting of cusp points:
152         k = find(d(K).*d(K+1)>0)
153         K(k) = K(k)+1
154
155         // We unfold:
156         // The level and the slope at the left edge are taken as references.
157         // The slope may be not of convenient sign. In this case, the
158         //  unfolded profile will be upside down
159         // One section over two must be kept as is, others must be flipped
160
161         nCusp = length(K)   // number of detected cusp points
162         if nCusp>0 then
163             flip = %t
164             // Main loop over cusp points:
165             for p = 1:nCusp
166                 if flip then
167                     iLeft = K(p)+1
168                     if p<nCusp then
169                         iRight = K(p+1)
170                     else
171                         iRight = nC
172                     end
173                     tmp = iLeft:iRight
174                     d(tmp) = -d(tmp)        // the section is flipped
175                     d(iLeft) = (d(iLeft-1)+d(iLeft+1))/2 // the reconnexion is smoothed
176                 end
177                 flip = ~flip    // one section / 2 is left as is
178             end
179         end
180
181         // trimming both edges =>  n-1 columns:
182         d = d(:, 2:$-1)
183
184         // The first point is taken as anchor:
185         r = [a(1)  d]
186     end
187
188     // Building the new profile from its restored derivative
189     retval = cumsum(r,"c")
190
191     // Post-processing
192     if transposed then
193         retval = retval.'
194     end
195 endfunction