* Bug #9819 fixed - Elementary_functions: unwrap 57/13357/5
Samuel Gougeon [Wed, 11 Dec 2013 09:12:56 +0000 (10:12 +0100)]
unwrap did not exist in Scilab.

Change-Id: I39a2dc18db577209359d2ca351b6c80a5a761377

17 files changed:
SEP/INDEX
SEP/SEP_126_unwrap.odt [new file with mode: 0644]
scilab/CHANGES_5.5.X
scilab/modules/elementary_functions/help/en_US/modulo.xml
scilab/modules/elementary_functions/help/en_US/trigonometry/atan.xml
scilab/modules/elementary_functions/help/en_US/unwrap.xml [new file with mode: 0644]
scilab/modules/elementary_functions/help/fr_FR/modulo.xml
scilab/modules/elementary_functions/help/ja_JP/trigonometry/atan.xml
scilab/modules/elementary_functions/help/pt_BR/trigonometry/atan.xml
scilab/modules/elementary_functions/help/ru_RU/trigonometry/atan.xml
scilab/modules/elementary_functions/macros/%_unwrap.sci [new file with mode: 0644]
scilab/modules/elementary_functions/macros/unwrap.sci [new file with mode: 0644]
scilab/modules/elementary_functions/tests/unit_tests/unwrap.dia.ref [new file with mode: 0644]
scilab/modules/elementary_functions/tests/unit_tests/unwrap.tst [new file with mode: 0644]
scilab/modules/helptools/etc/images_md5.txt
scilab/modules/helptools/images/unwrap_1.png [new file with mode: 0644]
scilab/modules/helptools/images/unwrap_2.png [new file with mode: 0644]

index 5ccb3de..6ea00fa 100644 (file)
--- a/SEP/INDEX
+++ b/SEP/INDEX
@@ -121,3 +121,4 @@ SEP #122: Customized mark draw
 SEP #123: Grid improvements
 SEP #124: New input argument for csvRead(), to ignore header comments.
 SEP #125: New input argument to matfile_open to open Matlab 7.3 files.
+SEP #126: New function unwrap
diff --git a/SEP/SEP_126_unwrap.odt b/SEP/SEP_126_unwrap.odt
new file mode 100644 (file)
index 0000000..3f937b9
Binary files /dev/null and b/SEP/SEP_126_unwrap.odt differ
index e5b913e..caf2739 100644 (file)
@@ -11,6 +11,7 @@ New Features
  - jcreatejar - Creates a Java archive (JAR) from a set of files / directories
  - ilib_build_jar - Builds Java packages from sources into a JAR file
  - ifftshift - inverse FFT shift
+ - unwrap - unwraps/unfolds a Y(x) profile
 
 * modulo() and pmodulo() now support integers & hypermatrices (See bug #13002).
 
@@ -156,6 +157,8 @@ Scilab Bug Fixes
 
 * Bug #9701 fixed - optim with "qn" option was failing for large problems.
 
+* Bug #9819 fixed - unwrap function did not exist in Scilab.
+
 * Bug #9840 fixed - Default G tolerance in lsqrsolve was too large.
 
 * Bug #10012 fixed - lmisolver and lmitool called with no input now produce errors.
index 48d49f0..f376978 100644 (file)
@@ -3,7 +3,7 @@
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2008 - INRIA
  * Copyright (C) 2013 - Samuel GOUGEON
- * 
+ *
  * This file must be used under the terms of the CeCILL.
  * This source file is licensed as described in the file COPYING, which
  * you should have received as part of this distribution.  The terms
@@ -31,8 +31,8 @@
             <varlistentry>
                 <term>m, n</term>
                 <listitem>
-                    <para>Scalar, vector, matrix or hypermatrix of encoded integers, reals or polynomials (Hypermatrix is not supported for polynomials). 
-                        <varname>m</varname> and <varname>n</varname> must have the same type. If they are of integer type, they may be of distinct 
+                    <para>Scalar, vector, matrix or hypermatrix of encoded integers, reals or polynomials (Hypermatrix is not supported for polynomials).
+                        <varname>m</varname> and <varname>n</varname> must have the same type. If they are of integer type, they may be of distinct
                         encoding length (for instance int8 and int16). If none of them is scalar, they must have the same sizes.
                     </para>
                 </listitem>
@@ -68,7 +68,7 @@
     </refsection>
     <refsection>
         <title>Examples</title>
-        <programlisting role="example"><![CDATA[ 
+        <programlisting role="example"><![CDATA[
 n = [1,2,10,15];
 m = [2,2,3,5];
 modulo(n,m)
@@ -112,6 +112,9 @@ pmodulo(%z^2+1, %z)
         <title>See also</title>
         <simplelist type="inline">
             <member>
+                <link linkend="unwrap">unwrap</link>
+            </member>
+            <member>
                 <link linkend="ieee">ieee</link>
             </member>
         </simplelist>
index 0a20cf5..53cfbf4 100644 (file)
@@ -66,7 +66,7 @@
     </refsection>
     <refsection>
         <title>Examples</title>
-        <programlisting role="example"><![CDATA[ 
+        <programlisting role="example"><![CDATA[
 // examples with the second form
 x=[1,%i,-1,%i]
 phase_x=atan(imag(x),real(x))
@@ -92,6 +92,9 @@ atan(-%i)
                 <link linkend="tan">tan</link>
             </member>
             <member>
+                <link linkend="unwrap">unwrap</link>
+            </member>
+            <member>
                 <link linkend="ieee">ieee</link>
             </member>
         </simplelist>
diff --git a/scilab/modules/elementary_functions/help/en_US/unwrap.xml b/scilab/modules/elementary_functions/help/en_US/unwrap.xml
new file mode 100644 (file)
index 0000000..6267e42
--- /dev/null
@@ -0,0 +1,300 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!--
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) 2013 - Samuel GOUGEON
+ *
+ * This file must be used under the terms of the CeCILL.
+ * This source file is licensed as described in the file COPYING,
+ * which you should have received as part of this distribution.
+ * The terms are also available at
+ * http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+ *
+ -->
+<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns5="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="unwrap" xml:lang="en">
+    <refnamediv>
+        <refname>unwrap</refname>
+        <refpurpose>unwrap a Y(x) profile or a Z(x,y) surface. Unfold a Y(x) profile</refpurpose>
+    </refnamediv>
+    <refsynopsisdiv>
+        <title>Calling Sequences</title>
+        <synopsis>unwrap() // runs some examples
+            [U, breakPoints] = unwrap(Y)
+            [U, breakPoints] = unwrap(Y, z_jump)
+            [U, cuspPoints]  = unwrap(Y, "unfold")
+            U = unwrap(Z)
+            U = unwrap(Z, z_jump)
+            U = unwrap(Z, z_jump, dir)
+        </synopsis>
+    </refsynopsisdiv>
+    <refsection>
+        <title>Arguments</title>
+        <variablelist>
+            <varlistentry>
+                <term>Y</term>
+                <listitem>
+                    <para>Vector of real numbers: the profile to unwrap or unfold. Implicit abscissae X are assumed to be equispaced.</para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>Z</term>
+                <listitem>
+                    <para>Matrix of real numbers: the surface to unwrap. Implicit abscissae (X,Y) are assumed to be cartesian and equispaced (constant steps may be different along X versus along Y).</para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>z_jump</term>
+                <listitem>
+                    <para>
+                        Scalar real positive number used in unwrapping mode: the jump's height applied at breakpoints, performing the unwrapping. Only its absolute value is considered. The jump actually applied has the sign of the slope on both sides of each breakpoint. The default value is <literal>z_jump = 2*%pi</literal>. The special value <literal>z_jump = 0</literal> applies jumps equal to the average slope around each breakpoint, restoring a continuous slope over the whole profile or surface.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>dir</term>
+                <listitem>
+                    <para>"c" | "r" | "" (default): direction along which unwrapping is performed. "c" unwraps along columns, "r" unwraps along rows, "" unwraps in both directions.</para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>"unfold"</term>
+                <listitem>
+                    <para>Provide this switch to unfold the given curve if it is folded, instead of unwrapping it.</para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>U</term>
+                <listitem>
+                    <para>
+                        Unwrapped profile or surface, or unfolded profile. <varname>U</varname> has the same sizes as <varname>Y</varname> or <varname>Z</varname>.
+                    </para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>breakPoints, cuspPoints</term>
+                <listitem>
+                    <para>
+                        Vector of indices of points in <varname>Y</varname> where wrapping or folding has been detected and processed.
+                    </para>
+                </listitem>
+            </varlistentry>
+        </variablelist>
+    </refsection>
+    <refsection>
+        <title>Description</title>
+        <para>UNWRAPPING</para>
+        <para>
+            <function>unwrap()</function> will be useful to process profiles or even surfaces wrapped for instance by a periodic and monotonic function such as
+            <literal>Y = modulo(X,w)</literal> or <literal>Y = atan(tan(X))</literal>. It aims to somewhat invert these functions, recovering the input <literal>X</literal> over it full range instead of the limited <literal>w</literal> or <literal>[-%pi/2, %pi/2]</literal> one.
+        </para>
+        <para>A breakpoint of a wrapped profile is detected as a point where slopes on both neighbouring sides of the point are almost equal but much smaller (in absolute value) from and opposite to the slope at the considered point: at the point, there is a jump breaking and opposed to the neighbouring slope.
+        </para>
+        <para>This detection strategy avoids considering any particular level as a wrapping one. It allows to process wrapped profiles to which a constant (or even a trend) has been added afterwards.
+        </para>
+        <para>
+            Unwrapping consists in reducing every detected jump and somewhat restoring a continuous slope (initially assumed to be so). At each point, it is performed by applying a Y-shift on a whole side of the profile, with respect to the other. The Y-shift may be the same for all breakpoints, as provided by the user. If the user specifies a null Y-shift, <function>unwrap()</function> applies a jump equal to the average neighbouring slope, depending on each breakpoint.
+        </para>
+        <warning>An unwrapped profile is always defined but for a constant.</warning>
+        <para>
+            Unless <varname>dir</varname> is used, <function>unwrap()</function> applied to a surface unwraps its first column. Each point of this one is then taken as reference level for unwrapping each line starting with it.
+        </para>
+        <para> </para>
+        <para>UNFOLDING</para>
+        <para>
+            If the <varname>"unfold"</varname> keyword is used and a profile -- not a surface -- is provided, the profile is assumed to be folded instead of being wrapped.
+        </para>
+        <para>At a folding point -- or "cusp point" --, the profile is continuous, but its slope is broken: the slope has almost the same absolute value on both sides of the point, but is somewhat opposed from one side to the other.
+        </para>
+        <para>
+            Folding may occur for instance when a non-monotonic periodic function and its inverse are applied to a profile X, like with <varname>Y= acos(cos(X))</varname>. Recovering X from Y is quite more difficult than if it was wrapped. Indeed, some ambiguous situations may exist, like when the profile is tangentially folded on one of its quasi-horizontal sections (if any).
+        </para>
+        <para>When a cusp point is detected, a) one side of the profile starting from it is opposed (upside down), and b) the continuity of the profile and of its slope are preserved and retrieved at the considered point (this may need adding a small jump by the local slope).</para>
+        <warning>The slope on the left edge of the input profile is used as starting reference. The unfolded profile may be upside down with respect to the original true one. In addition, as for unwrapping, it is defined but for a constant.
+        </warning>
+        <para>Known limitations: presently, folded surfaces can't be processed.</para>
+    </refsection>
+    <refsection>
+        <title>Examples</title>
+        <para>Unwrapping or unfolding 1D profiles:</para>
+        <programlisting role="example"><![CDATA[
+// 1D EXAMPLES
+// -----------
+f = scf();
+f.figure_size = [800 1000];
+f.figure_position(2) = 0;
+f.figure_name = _("unwrap() & ""unfold"": 1D examples ");
+ax = gda();
+ax.y_label.font_size=2;
+drawlater()
+
+// Original 1D profile
+t = linspace(-4,4.2,800);
+alpha = t.^2 + t -1;
+subplot(5,2,1)
+titlepage("unwrap(): unwrap | unfold")
+subplot(5,2,2)
+plot(t,alpha)
+t2 = "$\text{Original profile: } \alpha=t^2+t-1$";
+ax = gca();
+ax.tight_limits = "on";
+yT = max(alpha) - strange(alpha)*0.1;
+xstring(0,yT,t2)
+e = gce();
+e.text_box_mode = "centered";
+e.font_size = 2;
+
+// Loop over cases
+for i=1:4
+    subplot(5,2,2*i+1)
+    if i==1 then
+        // Wrapping by atan(tan())
+        ralpha = atan(tan(alpha));            // Raw recovered alpha [pi]
+        ylabel("$atan(tan(\alpha))$")
+        [u, K] = unwrap(ralpha, %pi);         // arctan
+        t2 = "$\text{unwrap(Y, \%pi)}$";
+    elseif i==2
+        // Wrapping by modulo() + Y-shift
+        c = (rand(1,1)-0.5)*4;
+        ralpha = pmodulo(alpha, 5) + c;
+        ylabel("$modulo(\alpha,\ 5)"+sprintf("%+5.2f",c)+"$")
+        [u, K] = unwrap(ralpha, 0);
+        t2 = "$\text{unwrap(Y, 0)}$";
+    elseif i==3
+        // Folding by asin(sin()) + Y-shift
+        ralpha = 1+asin(sin(alpha));          // Raw recovered alpha [2.pi]
+        ylabel("$1+asin(sin(\alpha))$")
+        [u, K] = unwrap(ralpha, "unfold");
+        t2 = "$\text{unwrap(Y,""unfold"")}$";
+    else
+        // Folding by acos(cos()) + a trend
+        ralpha = 1+alpha/10+acos(cos(alpha)); // Raw recovered alpha [2.pi]
+        ylabel("$1+\frac{\alpha}{10}+acos(cos(\alpha))$")
+        [u, K] = unwrap(ralpha, "unfold");
+        t2 = "$\text{unwrap(Y,""unfold"")}$";
+    end
+    // Plotting the profile to be processed
+    plot(t, ralpha)
+    // Staring the breakpoints or the cusp points on the curve:
+    if K~=[] then
+        plot(t(K), ralpha(K),"*")
+    end
+    // Plotting the processed (unwrapped/unfolded) profile:
+    subplot(5,2,2*i+2)
+    plot(t,u)
+    ax = gca();
+    ax.tight_limits = "on";
+    // Adding a legend:
+    yT = max(u) - strange(u)*0.2;
+    xstring(0,yT,t2)
+    e = gce();
+    e.text_box_mode = "centered";
+    e.font_size = 2;
+end
+sda();
+drawnow()
+        ]]></programlisting>
+        <scilab:image>
+            %_unwrap()
+        </scilab:image>
+        <para>Unwrapping 2D surfaces:</para>
+        <programlisting role="example"><![CDATA[
+// 2D EXAMPLES
+// -----------
+ax = gda();
+ax.title.font_size = 2;
+f = scf();
+f.color_map = hotcolormap(100);
+f.figure_size = [475 1050];
+f.figure_position(2) = 0;
+f.figure_name = _("unwrap(): 2D examples");
+drawlater()
+
+nx = 300;
+ny = 400;
+rmax = 8.8;
+x = linspace(-rmax/2, rmax/2, nx)-1;
+y = linspace(-rmax/2, rmax/2, ny)+1;
+[X, Y] = meshgrid(x,y);
+for ex=0:1   // examples
+    // Original surface
+        // Generating the surface
+    if ex==0 then
+        z = X.^2 + Y.^2;
+    else
+        z = sqrt(0.3+sinc(sqrt(z)*3))*17-7;
+    end
+    // PLotting it in 3D
+    subplot(4,2,1+ex)
+    surf(x, y, z)
+    title("Original profile Z")
+    e = gce();
+    e.thickness = 0; // removes the mesh
+    e.parent.tight_limits = "on";
+
+    // Wrapped surface (flat)
+    m = 2.8;
+    zw = pmodulo(z, m); // wraps it
+    subplot(4,2,3+ex)
+    grayplot(x, y, zw')
+    title(msprintf("Zw = pmodulo(Z, %g)  (flat)",m))
+    ax0 = gca();
+    ax0.tight_limits = "on";
+
+    // Unwrapped surfaces (flat):
+    // in both directions:
+    u = unwrap(zw, 0);
+    subplot(4,2,5+ex)
+    grayplot(x, y, u')
+    title(msprintf("unwrap(Zw, %g)  (flat)", 0))
+    ax3 = gca();
+    ax3.tight_limits = "on";
+
+    if ex==0 then
+        direc = "r";
+    else
+        direc = "c";
+    end
+    // Along a single direction:
+    u = unwrap(zw, m, direc);
+    subplot(4,2,7+ex)
+    grayplot(x, y, u')
+    title(msprintf("unwrap(Zw, %g, ""%s"")  (flat)",m,direc))
+    ax1 = gca();
+    ax1.tight_limits = "on";
+end
+sda();
+drawnow()
+        ]]></programlisting>
+        <scilab:image>
+            %_unwrap("2D")
+        </scilab:image>
+    </refsection>
+    <refsection role="see also">
+        <title>See Also</title>
+        <simplelist type="inline">
+            <member>
+                <link linkend="atan">atan</link>
+            </member>
+            <member>
+                <link linkend="modulo">modulo</link>
+            </member>
+        </simplelist>
+    </refsection>
+    <refsection role="bibliography">
+        <title>Bibliography</title>
+        <para>
+            <ulink url="http://siptoolbox.sourceforge.net/doc/sip-0.7.0-reference/unwrapl.html">SIP toolbox: unwrap linear</ulink>
+        </para>
+        <para>
+            <ulink url="http://siptoolbox.sourceforge.net/doc/sip-0.7.0-reference/unwrapp.html">SIP toolbox: unwrap along path</ulink>
+        </para>
+    </refsection>
+    <refsection role="history tag">
+        <title>History</title>
+        <revhistory>
+            <revision>
+                <revnumber>5.5.0</revnumber>
+                <revdescription>unwrap() function introduced</revdescription>
+            </revision>
+        </revhistory>
+    </refsection>
+</refentry>
index fd25312..23313b3 100644 (file)
@@ -3,7 +3,7 @@
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2008 - INRIA
  * Copyright (C) 2013 - Samuel GOUGEON
- * 
+ *
  * This file must be used under the terms of the CeCILL.
  * This source file is licensed as described in the file COPYING, which
  * you should have received as part of this distribution.  The terms
@@ -31,7 +31,7 @@
             <varlistentry>
                 <term>m, n</term>
                 <listitem>
-                    <para>Scalaire, vecteur, matrice ou hypermatrice d'entiers encodés, de décimaux réels ou de polynômes (les hypermatrices de polynômes ne sont pas admises). 
+                    <para>Scalaire, vecteur, matrice ou hypermatrice d'entiers encodés, de décimaux réels ou de polynômes (les hypermatrices de polynômes ne sont pas admises).
                         <varname>m</varname> et <varname>n</varname> doivent être de même type. S'ils sont de types entiers, il peuvent être d'encodages distincts (par exemple int8 et int16) Si aucune des deux n'est scalaire, elles doivent avoir les mêmes dimensions.
                     </para>
                 </listitem>
@@ -66,7 +66,7 @@
     </refsection>
     <refsection>
         <title>Exemples</title>
-        <programlisting role="example"><![CDATA[ 
+        <programlisting role="example"><![CDATA[
 n = [1,2,10,15];
 m = [2,2,3,5];
 modulo(n,m)
@@ -109,6 +109,9 @@ pmodulo(%z^2+1, %z)
         <title>Voir aussi</title>
         <simplelist type="inline">
             <member>
+                <link linkend="unwrap">unwrap</link>
+            </member>
+            <member>
                 <link linkend="ieee">ieee</link>
             </member>
         </simplelist>
index 1c08843..49d081d 100644 (file)
@@ -70,7 +70,7 @@
     </refsection>
     <refsection>
         <title>例</title>
-        <programlisting role="example"><![CDATA[ 
+        <programlisting role="example"><![CDATA[
 // 二番目の形式の例
 x=[1,%i,-1,%i]
 phasex=atan(imag(x),real(x))
@@ -96,6 +96,9 @@ atan(-%i)
                 <link linkend="tan">tan</link>
             </member>
             <member>
+                <link linkend="unwrap">unwrap</link>
+            </member>
+            <member>
                 <link linkend="ieee">ieee</link>
             </member>
         </simplelist>
index 0e5a258..66cd790 100644 (file)
@@ -69,7 +69,7 @@
     </refsection>
     <refsection>
         <title>Exemplos</title>
-        <programlisting role="example"><![CDATA[ 
+        <programlisting role="example"><![CDATA[
 // exemplos com a segunda forma
 x=[1,%i,-1,%i]
 phasex=atan(imag(x),real(x))
@@ -95,6 +95,9 @@ atan(-%i)
                 <link linkend="tan">tan</link>
             </member>
             <member>
+                <link linkend="unwrap">unwrap</link>
+            </member>
+            <member>
                 <link linkend="ieee">ieee</link>
             </member>
         </simplelist>
index 690688b..035abef 100644 (file)
@@ -52,7 +52,7 @@
             Диапазон <code>atan(y, x)</code> равен <latex>(-\pi, \pi]</latex>.
         </para>
         <para>
-            Для вещественных аргументов обе формы дают идентичные значения, если 
+            Для вещественных аргументов обе формы дают идентичные значения, если
             <literal>x&gt;0</literal>.
         </para>
         <para>
     </refsection>
     <refsection>
         <title>Примеры</title>
-        <programlisting role="example"><![CDATA[ 
+        <programlisting role="example"><![CDATA[
 // примеры со второй формой
 x=[1,%i,-1,%i]
 phase_x = atan(imag(x),real(x))
 atan(0,-1)
 atan(-%eps,-1)
+
 // переходы
 atan(-%eps + 2*%i)
 atan(+%eps + 2*%i)
 atan(-%eps - 2*%i)
 atan(+%eps - 2*%i)
+
 // значения в точках перехода
 ieee(2)
 atan(%i)
@@ -89,6 +89,9 @@ atan(-%i)
                 <link linkend="tan">tan</link>
             </member>
             <member>
+                <link linkend="unwrap">unwrap</link>
+            </member>
+            <member>
                 <link linkend="ieee">ieee</link>
             </member>
         </simplelist>
diff --git a/scilab/modules/elementary_functions/macros/%_unwrap.sci b/scilab/modules/elementary_functions/macros/%_unwrap.sci
new file mode 100644 (file)
index 0000000..bec292e
--- /dev/null
@@ -0,0 +1,155 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) - 2013 - Samuel GOUGEON
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function %_unwrap(typexample)
+    if argn(2)==0 then
+        typexample = ""
+    end
+    if convstr(typexample)~="2d" then
+        // 1D EXAMPLES
+        // -----------
+        f = scf();
+        f.figure_size = [800 1000];
+        f.figure_position(2) = 0;
+        f.figure_name = _("unwrap() & ""unfold"": 1D examples ");
+        ax = gda();
+        ax.y_label.font_size=2;
+        drawlater()
+
+        // Original 1D profile
+        t = linspace(-4,4.2,800);
+        alpha = t.^2 + t -1;
+        subplot(5,2,1)
+        titlepage("unwrap(): unwrap | unfold")
+        subplot(5,2,2)
+        plot(t,alpha)
+        t2 = "$\text{Original profile: } \alpha=t^2+t-1$";
+        ax = gca();
+        ax.tight_limits = "on";
+        yT = max(alpha) - strange(alpha)*0.1;
+        xstring(0,yT,t2)
+        e = gce();
+        e.text_box_mode = "centered";
+        e.font_size = 2;
+
+        // Loop over cases
+        for i=1:4
+            subplot(5,2,2*i+1)
+            if i==1 then
+                // Wrapping by atan(tan())
+                ralpha = atan(tan(alpha));            // raw recovered alpha [pi]
+                ylabel("$\mbox{Y= }atan(tan(\alpha))$")
+                [u, K] = unwrap(ralpha, %pi);         // arctan
+                t2 = "$\text{unwrap(Y, \%pi)}$";
+            elseif i==2
+                // Wrapping by modulo() + Y-shift
+                c = (rand(1,1)-0.5)*4;
+                ralpha = pmodulo(alpha, 5) + c;
+                ylabel("$\mbox{Y= }modulo(\alpha,\ 5)"+sprintf("%+5.2f",c)+"$")
+                [u, K] = unwrap(ralpha, 0);
+                t2 = "$\text{unwrap(Y, 0)}$";
+            elseif i==3
+                // Folding by asin(sin()) + Y-shift
+                ralpha = 1+asin(sin(alpha));          // raw recovered alpha [2.pi]
+                ylabel("$\mbox{Y= }1+asin(sin(\alpha))$")
+                [u, K] = unwrap(ralpha, "unfold");
+                t2 = "$\text{unwrap(Y,""unfold"")}$";
+            else
+                // Folding by acos(cos()) + trend
+                ralpha = 1+alpha/10+acos(cos(alpha)); // raw recovered alpha [2.pi]
+                ylabel("$\mbox{Y= }1+\frac{\alpha}{10}+acos(cos(\alpha))$")
+                [u, K] = unwrap(ralpha, "unfold");
+                t2 = "$\text{unwrap(Y,""unfold"")}$";
+            end
+            // Plotting the profile to be processed
+            plot(t, ralpha)
+            // Staring the breakpoints or the cusp points on the curve:
+            if K~=[]
+                plot(t(K), ralpha(K),"*")
+            end
+            // Plotting the processed (unwrapped/unfolded) profile:
+            subplot(5,2,2*i+2)
+            plot(t,u)
+            ax = gca();
+            ax.tight_limits = "on";
+            // Adding a legend:
+            yT = max(u) - strange(u)*0.2;
+            xstring(0,yT,t2)
+            e = gce();
+            e.text_box_mode = "centered";
+            e.font_size = 2;
+        end
+        sda();
+        drawnow()
+
+    else
+        // 2D EXAMPLES
+        // -----------
+        ax = gda();
+        ax.title.font_size = 2;
+        f = scf();
+        f.color_map = hotcolormap(100);
+        f.figure_size = [475 1050];
+        f.figure_position(2) = 0;
+        f.figure_name = _("unwrap(): 2D examples");
+        drawlater()
+
+        nx = 300;
+        ny = 400;
+        rmax = 8.8;
+        x = linspace(-rmax/2, rmax/2, nx)-1;
+        y = linspace(-rmax/2, rmax/2, ny)+1;
+        [X, Y] = meshgrid(x,y);
+        for ex=0:1   // examples
+            // Original surface
+            // Generating the surface
+            if ex==0 then
+                z = X.^2 + Y.^2;
+            else
+                z = sqrt(0.3+sinc(sqrt(z)*3))*17-7;
+            end
+            // PLotting it in 3D
+            subplot(4,2,1+ex)
+            surf(x, y, z)
+            title("Original profile Z")
+            e = gce();
+            e.thickness = 0; // removes the mesh
+            e.parent.tight_limits = "on";
+
+            // WRAPPED surface (flat)
+            m = 2.8;
+            zw = pmodulo(z, m); // wraps it
+            subplot(4,2,3+ex)
+            grayplot(x, y, zw')
+            title(msprintf("Zw = pmodulo(Z, %g)  (flat)",m))
+            ax0 = gca();
+            ax0.tight_limits = "on";
+
+            // UNWRAPPED surfaces (flat):
+            // in both directions:
+            u = unwrap(zw, 0);
+            subplot(4,2,5+ex)
+            grayplot(x, y, u')
+            title(msprintf("unwrap(Zw, %g)  (flat)", 0))
+            ax3 = gca();
+            ax3.tight_limits = "on";
+
+            if ex==0, direc="r", else direc="c", end
+            // along a single direction:
+            u = unwrap(zw, m, direc);
+            subplot(4,2,7+ex)
+            grayplot(x, y, u')
+            title(msprintf("unwrap(Zw, %g, ""%s"")  (flat)",m,direc))
+            ax1 = gca();
+            ax1.tight_limits = "on";
+        end
+        sda();
+        drawnow()
+    end
+endfunction
diff --git a/scilab/modules/elementary_functions/macros/unwrap.sci b/scilab/modules/elementary_functions/macros/unwrap.sci
new file mode 100644 (file)
index 0000000..40a1ae9
--- /dev/null
@@ -0,0 +1,190 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) - 2013 - Samuel GOUGEON
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [retval, K] = unwrap(a, varargin)
+    // [unwr, cuspPoints] = unwrap(a)        // row or column vector or matrix
+    // [unwr, cuspPoints] = unwrap(a, jump)  // default jump = 2*%pi
+    // jump=0 => = average local slope
+    // [unwr, cuspPoints] = unwrap(a, jump, direc)
+    //          direc="c": unwrap/unfold along columns
+    //          direc="r": along rows
+    //          else: along both directions (default)
+    // [unwr, cuspPoints] = unwrap(a, "unfold")
+
+    // Initializations
+    // --------------
+    [lhs,rhs] = argn(0);
+    transposed = %f
+    unfold = %f     // unfold (1D) instead of unwrap (1D|2D)
+    jump = 2*%pi
+    direc  = ""
+    retval = []
+    K = []
+
+    // EXAMPLES
+    // --------
+    if rhs==0 then
+        %_unwrap()      // display 1D examples
+        halt(_("Press return to display 2D examples"))
+        %_unwrap("2D")  // display 2D examples
+        return
+    end
+
+    // =================== CHECKING PARAMETERS ====================
+    if type(a)~=1 then
+        msg = _("%s: Wrong type for input argument #%d: Real expected.\n")
+        error(msprintf(msg, "unwrap",1))
+    end
+    if rhs>1
+        if typeof(varargin(1))=="string" then
+            if varargin(1)=="unfold" then
+                unfold = %t
+            else
+                msg = _("%s: Wrong value for input argument #%d: Must be in the set  {%s}.\n")
+                error(msprintf(msg, "unwrap",2,"""unfold"""))
+            end
+        else
+            jump = varargin(1)
+            if type(jump)~=1 then
+                msg = _("%s: Wrong type for input argument #%d: Real expected.\n")
+                error(msprintf(msg, "unwrap",2))
+            else
+                jump = abs(jump(1))
+            end
+        end
+        if ~unfold & rhs>=3 then
+            direc = varargin(2)
+            if typeof(direc)~="string" then
+                msg = _("%s: Wrong type for input argument #%d: String expected.\n")
+                error(msprintf(msg, "unwrap", 3))
+            else
+                direc = direc(1)
+                if and(direc~=["r" "c"]) then
+                    msg = _("%s: Wrong value for input argument #%d: Must be in the set  {%s}.\n")
+                    error(msprintf(msg, "unwrap", 3, """r"",""c"""))
+                end
+            end
+        end
+    end
+
+    if direc=="c" | size(a,2)==1
+        a = a.'
+        transposed = %t
+    end
+
+    [m,n] = size(a)
+    if (n < 4)
+        if transposed then
+            msg = _("%s: not enough tall input argument #%d: must have at least 4 rows\n")
+        else
+            msg = _("%s: too narrow input argument #%d: must have at least 4 columns\n")
+        end
+        error(msprintf(msg, "unwrap",1))
+    end
+    if unfold & m>1 then
+        msg = _("%s: Wrong size for input argument #%d: Vector expected.\n")
+        error(msprintf(msg, "unwrap",1))
+    end
+
+    // ============================ PROCESSING ==============================
+
+    // Columns #1 and #$ are duplicated:
+    rae = [a(:,1)  a  a(:,$)];       // => n+2 columns
+    // Derivative along rows (assuming equally spaced x values)
+    d = rae(:,2:$) - rae(:, 1:$-1);     // n+1 columns
+    clear rae
+
+    // ---------------------------  U N W R A P  ---------------------------
+    if ~unfold then
+        // jump (may be)  (n-1 columns)
+        ju = d(:, 2:$-1)
+        // average Local derivative (before/after jump) (n-1 columns) :
+        avL = (d(:, 1:$-2)+d(:, 3:$))/2
+        // where
+        wh = abs(ju)>4*abs(avL) & d(:, 1:$-2).*d(:, 3:$)>0
+        wh = abs(ju)>5*abs(avL)
+        K = find(wh)
+        d = d(:, 2:$-1)                  // n-1 columns
+        if jump~=0 then
+            d(K) = ju(K) - sign(ju(K)).*jump
+            // Cleaning wrongly inserted jumps :
+            k = find(abs(d(K)-mean(d))>5*stdev(d))
+            d(K(k)) = d(K(k)) - sign(d(K(k))).*jump
+        else
+            d(K) = avL(K)
+        end
+
+        if and(size(a)~=1) & direc==""
+            // 2D case: levels all rows according to the unwrapped 1st column
+            p = unwrap(a(:,1),jump)
+        else
+            p  = a(:,1)
+        end
+        r = [p d]
+
+        // ---------------------------  U N F O L D  ---------------------------
+    else
+        // Local radius of curvature: calculation
+        Rc = ((1+d(:,2:$).^2).^1.5) ./ abs(d(:,2:$)-d(:,1:$-1))
+        K = find(Rc < mean(Rc)/30) // criterium to detect cusp points
+        [nR,nC] = size(Rc)
+
+        // trimming edges (due to duplication of 1st & last columns):
+        tmp = find(K==1 | K==nC)
+        K(tmp) = []
+
+        // trimming duplicates (if any)
+        tmp = find(K(2:$)-K(1:$-1)==1)
+        K(tmp+1) = []
+
+        // fine targetting of cusp points:
+        k = find(d(K).*d(K+1)>0)
+        K(k) = K(k)+1
+
+        // We unfold:
+        // The level and the slope at the left edge are taken as references.
+        // The slope may be not of convenient sign. In this case, the
+        //  unfolded profile will be upside down
+        // One section over two must be kept as is, others must be flipped
+
+        nCusp = length(K)   // number of detected cusp points
+        if nCusp>0 then
+            flip = %t
+            // Main loop over cusp points:
+            for p = 1:nCusp
+                if flip then
+                    iLeft = K(p)+1
+                    if p<nCusp then
+                        iRight = K(p+1)
+                    else
+                        iRight = nC
+                    end
+                    tmp = iLeft:iRight
+                    d(tmp) = -d(tmp)        // the section is flipped
+                    d(iLeft) = (d(iLeft-1)+d(iLeft+1))/2 // the reconnexion is smoothed
+                end
+                flip = ~flip    // one section / 2 is left as is
+            end
+        end
+
+        // trimming both edges =>  n-1 columns:
+        d = d(:, 2:$-1)
+
+        // The first point is taken as anchor:
+        r = [a(1)  d]
+    end
+
+    // Building the new profile from its restored derivative
+    retval = cumsum(r,"c")
+
+    // Post-processing
+    if transposed then
+        retval = retval.'
+    end
+endfunction
diff --git a/scilab/modules/elementary_functions/tests/unit_tests/unwrap.dia.ref b/scilab/modules/elementary_functions/tests/unit_tests/unwrap.dia.ref
new file mode 100644 (file)
index 0000000..2fa5345
--- /dev/null
@@ -0,0 +1,77 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Samuel Gougeon
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+// Example #1 (1D wrapped)
+t = linspace(-4, 4.2, 800);
+Y = t.^2 + t -1;
+Yw = atan(tan(Y));  // Raw wrapped recovered Y [pi]
+Yunw = unwrap(Yw, %pi);
+assert_checkequal(size(Yunw), size(Y));
+assert_checkequal(max(Yunw), 8.273629385640850486539);
+assert_checkequal(min(Yunw), -13.816370473381940797708);
+// Example #2 (1D wrapped)
+rand("seed", 0);
+c = (rand(1, 1)-0.5)*4;
+Yw = pmodulo(Y, 5) + c;  // unwrap() is unsensitive to a constant and even to a trend
+[Yunw, K] = unwrap(Yw, 0);
+assert_checkequal(size(Yunw), size(Y));
+assert_checkequal(max(Yunw), 9.6852994618564913764658);
+assert_checkequal(min(Yunw), -12.404700397166287473283);
+assert_checkequal(K, [15 98 233 450 585 668 734 791]);
+// Example #3 (1D folded)
+Yf = 1 + asin(sin(Y));  // Raw folded recovered Y [2.pi], + constant Y shift
+Yunf = unwrap(Yf, "unfold");
+assert_checkequal(size(Yunf), size(Y));
+assert_checkalmostequal(max(Yunf), 9.4279967527199310950436, [], 1e-14);
+assert_checkalmostequal(min(Yunf), -12.436144421661456505035, [], 1e-14);
+// Example #4 (1D folded)
+Yf = 1 + Y/10 + acos(cos(Y)); // Raw recovered Y [2.pi], + linear trend
+[Yunf, K] = unwrap(Yf, "unfold");
+assert_checkequal(size(Yunf), size(Y));
+assert_checkalmostequal(max(Yunf), 15.801371238872524926933, [], 1e-14);
+assert_checkalmostequal(min(Yunf), -6.1759877188984493301405, [], 1e-14);
+assert_checkequal(K, [24 75 138 233 451 546 609 660 704 743 779]);
+//Unwrapping 2D surfaces:
+nx = 300;
+ny = 400;
+rmax = 8.8;
+x = linspace(-rmax/2, rmax/2, nx)-1;
+y = linspace(-rmax/2, rmax/2, ny)+1;
+[X, Y] = meshgrid(x, y);
+// Example #5 : 2D wrapped
+Z = X.^2 + Y.^2;        // Paraboloïd
+m = 2.8;
+Zw = pmodulo(Z, m);     // >raps it
+Zunw  = unwrap(Zw, 0);  // Unwraps it in both directions, with free jumps
+assert_checkequal(size(Zunw), size(Z));
+assert_checkequal(max(Zunw), 19.120000000000032969183);
+assert_checkequal(min(Zunw), -39.199790375290433531);
+ZunwR = unwrap(Zw, 0, "r"); // Unwraps it along rows
+assert_checkequal(size(ZunwR), size(Z));
+assert_checkequal(max(ZunwR), 2.78738299382542109583);
+assert_checkequal(min(ZunwR), -28.995103166214036605197);
+ZunwC = unwrap(Zw, 0, "c"); // Unwraps it along columns
+assert_checkequal(size(ZunwC), size(Z));
+assert_checkequal(max(ZunwC), 20.398076084160148724322);
+assert_checkequal(min(ZunwC), -11.462196183827872530969);
+// Example #6 : 2D wrapped
+Z = sqrt(0.3+sinc(sqrt(Z)*3))*17-7; // Defines the surface
+m = 2.8;                     // Wrapping step (jump's amplitude)
+Zw = pmodulo(Z, m);          // Wraps it
+Zunw  = unwrap(Zw, m);       // Unwraps it in both directions, with given jump
+assert_checkequal(size(Zunw), size(Z));
+assert_checkequal(max(Zunw), 12.380638179897072603808);
+assert_checkequal(min(Zunw), -2.109245308884592162713);
+ZunwR = unwrap(Zw, m, "r"); // Unwraps it along rows
+assert_checkequal(size(ZunwR), size(Z));
+assert_checkequal(max(ZunwR), 12.380638179897072603808);
+assert_checkequal(min(ZunwR), -2.109245308884592162713);
+ZunwC = unwrap(Zw, m, "c"); // Unwraps it along columns
+assert_checkequal(size(ZunwC), size(Z));
+assert_checkequal(max(ZunwC), 12.380638179897072603808);
+assert_checkequal(min(ZunwC), -2.1092453088845903863557);
diff --git a/scilab/modules/elementary_functions/tests/unit_tests/unwrap.tst b/scilab/modules/elementary_functions/tests/unit_tests/unwrap.tst
new file mode 100644 (file)
index 0000000..ee800e8
--- /dev/null
@@ -0,0 +1,84 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Samuel Gougeon
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+
+// Example #1 (1D wrapped)
+t = linspace(-4, 4.2, 800);
+Y = t.^2 + t -1;
+Yw = atan(tan(Y));  // Raw wrapped recovered Y [pi]
+Yunw = unwrap(Yw, %pi);
+assert_checkequal(size(Yunw), size(Y));
+assert_checkequal(max(Yunw), 8.273629385640850486539);
+assert_checkequal(min(Yunw), -13.816370473381940797708);
+
+// Example #2 (1D wrapped)
+rand("seed", 0);
+c = (rand(1, 1)-0.5)*4;
+Yw = pmodulo(Y, 5) + c;  // unwrap() is unsensitive to a constant and even to a trend
+[Yunw, K] = unwrap(Yw, 0);
+assert_checkequal(size(Yunw), size(Y));
+assert_checkequal(max(Yunw), 9.6852994618564913764658);
+assert_checkequal(min(Yunw), -12.404700397166287473283);
+assert_checkequal(K, [15 98 233 450 585 668 734 791]);
+
+// Example #3 (1D folded)
+Yf = 1 + asin(sin(Y));  // Raw folded recovered Y [2.pi], + constant Y shift
+Yunf = unwrap(Yf, "unfold");
+assert_checkequal(size(Yunf), size(Y));
+assert_checkalmostequal(max(Yunf), 9.4279967527199310950436, [], 1e-14);
+assert_checkalmostequal(min(Yunf), -12.436144421661456505035, [], 1e-14);
+
+// Example #4 (1D folded)
+Yf = 1 + Y/10 + acos(cos(Y)); // Raw recovered Y [2.pi], + linear trend
+[Yunf, K] = unwrap(Yf, "unfold");
+assert_checkequal(size(Yunf), size(Y));
+assert_checkalmostequal(max(Yunf), 15.801371238872524926933, [], 1e-14);
+assert_checkalmostequal(min(Yunf), -6.1759877188984493301405, [], 1e-14);
+assert_checkequal(K, [24 75 138 233 451 546 609 660 704 743 779]);
+
+//Unwrapping 2D surfaces:
+nx = 300;
+ny = 400;
+rmax = 8.8;
+x = linspace(-rmax/2, rmax/2, nx)-1;
+y = linspace(-rmax/2, rmax/2, ny)+1;
+[X, Y] = meshgrid(x, y);
+
+// Example #5 : 2D wrapped
+Z = X.^2 + Y.^2;        // Paraboloïd
+m = 2.8;
+Zw = pmodulo(Z, m);     // >raps it
+Zunw  = unwrap(Zw, 0);  // Unwraps it in both directions, with free jumps
+assert_checkequal(size(Zunw), size(Z));
+assert_checkequal(max(Zunw), 19.120000000000032969183);
+assert_checkequal(min(Zunw), -39.199790375290433531);
+ZunwR = unwrap(Zw, 0, "r"); // Unwraps it along rows
+assert_checkequal(size(ZunwR), size(Z));
+assert_checkequal(max(ZunwR), 2.78738299382542109583);
+assert_checkequal(min(ZunwR), -28.995103166214036605197);
+ZunwC = unwrap(Zw, 0, "c"); // Unwraps it along columns
+assert_checkequal(size(ZunwC), size(Z));
+assert_checkequal(max(ZunwC), 20.398076084160148724322);
+assert_checkequal(min(ZunwC), -11.462196183827872530969);
+
+// Example #6 : 2D wrapped
+Z = sqrt(0.3+sinc(sqrt(Z)*3))*17-7; // Defines the surface
+m = 2.8;                     // Wrapping step (jump's amplitude)
+Zw = pmodulo(Z, m);          // Wraps it
+Zunw  = unwrap(Zw, m);       // Unwraps it in both directions, with given jump
+assert_checkequal(size(Zunw), size(Z));
+assert_checkequal(max(Zunw), 12.380638179897072603808);
+assert_checkequal(min(Zunw), -2.109245308884592162713);
+ZunwR = unwrap(Zw, m, "r"); // Unwraps it along rows
+assert_checkequal(size(ZunwR), size(Z));
+assert_checkequal(max(ZunwR), 12.380638179897072603808);
+assert_checkequal(min(ZunwR), -2.109245308884592162713);
+ZunwC = unwrap(Zw, m, "c"); // Unwraps it along columns
+assert_checkequal(size(ZunwC), size(Z));
+assert_checkequal(max(ZunwC), 12.380638179897072603808);
+assert_checkequal(min(ZunwC), -2.1092453088845903863557);
index d16b5a0..c7f3fbb 100644 (file)
@@ -1186,6 +1186,8 @@ title_1.png=43bf4550c7c9b2de76ddd4eddb479b2b
 titlepage_1.png=917be99c9df0964672daa45c284de115
 titlepage_fr_FR_1.png=f76157d643854ea8716178823dd2680e
 trans_1.png=a5ddde2263a28ffaba21606bcfc2232
+unwrap_1.png=b55bd8834f6ed873c3c614d5b81a7484
+unwrap_2.png=ea92aa6f961b2fd3c5d82a18529aecbd
 wavread_1.png=7f4fa7fd3215e86c6beac1b99ed5c107
 whitecolormap_1.png=5c902a640b651c5ec7bcfdeba604e3c0
 wiener_1.png=7cf1c15abff7e9f0c2ca9bd45882a5c
diff --git a/scilab/modules/helptools/images/unwrap_1.png b/scilab/modules/helptools/images/unwrap_1.png
new file mode 100644 (file)
index 0000000..4c8071b
Binary files /dev/null and b/scilab/modules/helptools/images/unwrap_1.png differ
diff --git a/scilab/modules/helptools/images/unwrap_2.png b/scilab/modules/helptools/images/unwrap_2.png
new file mode 100644 (file)
index 0000000..e69de29