File: movslice.m

package info (click to toggle)
octave 9.4.0-1
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid, trixie
  • size: 144,300 kB
  • sloc: cpp: 332,784; ansic: 77,239; fortran: 20,963; objc: 9,396; sh: 8,213; yacc: 4,925; lex: 4,389; perl: 1,544; java: 1,366; awk: 1,259; makefile: 648; xml: 189
file content (178 lines) | stat: -rw-r--r-- 6,567 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
########################################################################
##
## Copyright (C) 2018-2024 The Octave Project Developers
##
## See the file COPYRIGHT.md in the top-level directory of this
## distribution or <https://octave.org/copyright/>.
##
## This file is part of Octave.
##
## Octave is free software: you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <https://www.gnu.org/licenses/>.
##
########################################################################

## -*- texinfo -*-
## @deftypefn  {} {@var{slcidx} =} movslice (@var{N}, @var{wlen})
## @deftypefnx {} {[@var{slcidx}, @var{C}, @var{Cpre}, @var{Cpost}, @var{win}] =} movslice (@dots{})
## Generate indices to slice a vector of length @var{N} into windows
## of length @var{wlen}.
##
## The input @var{N} must be a positive integer.
##
## The moving window length input @var{wlen} can either be a scalar not equal
## to 1 or a 2-element array of integers.  For scalar values, if odd the window
## is symmetric and includes @w{@code{(@var{wlen} - 1) / 2}}@ elements on
## either side of the central element.  If @var{wlen} is even the window is
## asymmetric and has @w{@code{@var{wlen}/2}}@ elements to the left of the
## central element and @w{@code{@var{wlen}/2 - 1}}@ elements to the right of
## the central element.  When @var{wlen} is a 2-element array,
## @w{@code{[@var{nb}, @var{na}]}}, the window includes @var{nb} elements to
## the left of the current element and @var{na} elements to the right of the
## current element.
##
## The output @var{slcidx} is an array of indices of the slices that fit fully
## within the vector, where each column is an individual slice as the window
## moves from left to right.  The slices have @var{wlen} elements for scalar
## @var{wlen}, or @w{@code{@var{nb} + @var{na} + 1}}@ elements for array valued
## @var{wlen}.
##
## Optional output @var{C} is an row vector of window center positions where
## the window stays fully within the vector.
##
## Optional outputs @var{Cpre} and @var{Cpost} contain the vector elements at
## the start and end of the vector, respectively, that result in the window
## extending beyond the ends of the vector.
##
## Optional output @var{win} is a column vector with the same number of rows
## as @var{slcidx} that contains the moving window defined as a center
## relative position stencil.
##
## @seealso{movfun}
## @end deftypefn

function [slcidx, C, Cpre, Cpost, win] = movslice (N, wlen)

  if (nargin != 2)
    print_usage ();
  endif

  ## Validate N
  if (! (isscalar (N) && isindex (N)))
    error ("movslice: N must be a positive integer");
  endif

  ## Validate window length
  if (! (isnumeric (wlen) && all (wlen >= 0) && fix (wlen) == wlen))
    error ("Octave:invalid-input-arg",
           "movslice: WLEN must be a scalar or 2-element array of integers >= 0");
  endif
  if (isscalar (wlen))
    ## Check for proper window length
    if (wlen == 1)
      error ("Octave:invalid-input-arg", "movslice: WLEN must be > 1");
    endif
  elseif (numel (wlen) == 2)
    ## wlen = [nb, na].  No further validation here.
  else
    error ("Octave:invalid-input-arg",
           "movslice: WLEN must be a scalar or 2-element array of integers >= 0");
  endif
  if (max (wlen) > N)
    error ("Octave:invalid-input-arg", ...
           "movslice: window length WLEN (%d) must be shorter than length along DIM (%d)", ...
           max (wlen), N);
  endif

  if (isscalar (wlen))
    if (mod (wlen, 2) == 1)
      ## Symmetric window
      nb = na = (wlen - 1) / 2;
      wlen = [nb, na];
    else
      ## Asymmetric window
      nb = wlen / 2;
      na = nb - 1;
      wlen = [nb, na];
    endif
  endif

  Cpre  = 1:wlen(1);              # centers that can't fit the pre-window
  Cnf   = N - wlen(2) + 1;        # first center that can't fit the post-window
  Cpost = Cnf:N;                  # centers that can't fit centered post-window
  C     = (wlen(1) + 1):(Cnf - 1);
  ## Convert C to minimum unsigned integer array large enough to hold indices.
  ## This can save significant memory in resulting slcidx array over using a
  ## double (8 bytes).
  if (N <= 255)
    C = uint8 (C);
  elseif (N <= 65535)
    C = uint16 (C);
  elseif (N <= 4294967295)
    C = uint32 (C);
  else
    C = uint64 (C);
  endif
  win   = (-wlen(1):wlen(2)).';
  slcidx = C + win;

endfunction

%!assert (double (movslice (10, 2)), [1:9; 2:10])
%!assert (double (movslice (10, 9)), [1:9; 2:10].')

%!test
%! [sl, c, cpre, cpost, win] = movslice (10, 4);
%! assert (double (sl), [1:7; 2:8; 3:9; 4:10]);
%! assert (double (c), 3:9);
%! assert (cpre, 1:2);
%! assert (cpost, 10);
%! assert (win, [-2:1:1].');

%!test
%! [sl, c, cpre, cpost, win] = movslice (10, 10);
%! assert (double (sl), [1:10].');
%! assert (double (c), 6);
%! assert (cpre, 1:5);
%! assert (cpost, 7:10);
%! assert (win, [-5:1:4].');

## Verify uint output.  Don't test uint64 due to excessive memory usage.
%!test
%! [sl, c] = movslice (10, 10);
%! assert (class (sl), "uint8");
%! assert (class (c), "uint8");
%! [sl, c] = movslice (1000, 1000);
%! assert (class (sl), "uint16");
%! assert (class (c), "uint16");
%! [sl, c] = movslice (100000, 100000);
%! assert (class (sl), "uint32");
%! assert (class (c), "uint32");

## Test input validation
%!error <Invalid call> movslice ()
%!error <Invalid call> movslice (1)
%!error <N must be a positive integer> movslice ([1 2], 1)
%!error <N must be a positive integer> movslice (0, 1)
%!error <WLEN must be .* array of integers> movslice (1, {1})
%!error <WLEN must be .* array of integers .= 0> movslice (1, -1)
%!error <WLEN must be .* array of integers> movslice (1, 1.5)
%!error <WLEN must be . 1> movslice (1, 1)
%!error <WLEN must be a scalar or 2-element array> movslice (1, [1, 2, 3])
%!error <WLEN \(3\) must be shorter than length along DIM \(1\)>
%! movslice (1, 3);
%!error <WLEN \(4\) must be shorter than length along DIM \(1\)>
%! movslice (1, [4, 1]);
%!error <WLEN \(5\) must be shorter than length along DIM \(1\)>
%! movslice (1, [1, 5]);