text:
function tet_mesh_tet_neighbors ( prefix )

%*****************************************************************************80
%
%% MAIN is the main program for TET_MESH_TET_NEIGHBORS.
%
%  Discussion:
%
%    TET_MESH_TET_NEIGHBORS manages the tet mesh neighbor calculation.
%
%    A tet mesh of order 4 or order 10 may be used.
%
%  Usage:
%
%    tet_mesh_tet_neighbors ( 'prefix' )
%
%    where
%
%    * 'prefix'_nodes.txt contains nodal coordinates (not used by this function);
%    * 'prefix'_elements.txt the element definitions;
%    * 'prefix'_element_neighbors.txt contains the element neigbhbors.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 February 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TET_MESH_TET_NEIGHBORS:\n' );
  fprintf ( 1, '  MATLAB version:\n' );
  fprintf ( 1, '  Read a tet mesh dataset of TETRA_NUM\n' );
  fprintf ( 1, '  tetrahedrons, u sing 4 or 10 nodes.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compute tet mesh neighbors and write this\n' );
  fprintf ( 1, '  information to a file.\n' );
%
%  Argument 1 is the common file prefix.
%
  if ( 1 <= nargin )

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'TET_MESH_TET_NEIGHBORS:\n' );
    prefix = input ( '  Please enter the filename prefix:' );

  end
%
%  Create the filenames.
%
  element_filename = strcat ( prefix, '_elements.txt' );
  neighbor_filename = strcat ( prefix, '_element_neighbors.txt' );
%
%  Read the tetra data.
%
  [ tetra_order, tetra_num ] = i4mat_header_read ( element_filename );

  if ( tetra_order ~= 4 && tetra_order ~= 10 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TET_MESH_TET_NEIGHBORS - Fatal error!\n' );
    fprintf ( 1, '  Data is not for a 4 or 10 node tet mesh.\n' );
    error ( 'TET_MESH_TET_NEIGHBORS - Fatal error!' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the header of "%s".\n', element_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Tetrahedron order = %d\n', tetra_order );
  fprintf ( 1, '  Number of tetras  = %d\n', tetra_num );

  tetra_node = i4mat_data_read ( element_filename, tetra_order, ...
    tetra_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the data in "%s".\n', element_filename );
%
%  Detect and correct 0-based indexing.
%
  tetra_node = mesh_base_one ( tetra_order, tetra_num, tetra_node );
%
%  Print first 5.
%
  i4mat_transpose_print_some ( tetra_order, tetra_num, ...
    tetra_node, 1, 1, tetra_order, 5, '  First 5 tetrahedrons:' );
%
%  Compute the neighbor information.
%
  tetra_neighbor = tet_mesh_neighbor_tets ( tetra_order, tetra_num, ...
    tetra_node );

  i4mat_transpose_print_some ( 4, tetra_num, ...
    tetra_neighbor, 1, 1, 4, 5, '  First 5 neighbors:' );
%
%  Write the neighbor information to a file.
%  Create the output file name from the input file name.
%
  i4mat_write ( neighbor_filename, 4, tetra_num, tetra_neighbor );

  fprintf ( 1, '  Wrote the tetrahedron neighbor information to "%s".\n', ...
    neighbor_filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TET_MESH_TET_NEIGHBORS:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function column_num = file_column_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_COLUMN_COUNT counts the columns in the first line of a file.
%
%  Discussion:
%
%    The file is assumed to be a simple text file.
%
%    Most lines of the file are presumed to consist of COLUMN_NUM words,
%    separated by spaces.  There may also be some blank lines, and some
%    comment lines, which have a "#" in column 1.
%
%    The routine tries to find the first non-comment non-blank line and
%    counts the number of words in that line.
%
%    If all lines are blanks or comments, it goes back and tries to analyze
%    a comment line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the file.
%
%    Output, integer COLUMN_NUM, the number of columns in the file.
%
  FALSE = 0;
  TRUE = 1;
%
%  Open the file.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_COLUMN_COUNT - Error!' );
  end
%
%  Read one line, but skip blank lines and comment lines.
%  Use FGETL so we drop the newline character!
%
  got_one = FALSE;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

    elseif ( line(1) == '#' )

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

  if ( got_one == FALSE )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Warning!\n' );
    fprintf ( 1, '  The file does not seem to contain any data.\n' );
    column_num = -1;
    return;
  end

  column_num = s_word_count ( line );

  return
end
function row_num = file_row_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_ROW_COUNT counts the number of row records in a file.
%
%  Discussion:
%
%    Each input line is a "RECORD".
%
%    The records are divided into three groups:
%
%    * BLANK LINES (nothing but blanks)
%    * COMMENT LINES (begin with a '#')
%    * DATA RECORDS (anything else)
%
%    The value returned by the function is the number of data records.
%
%    By the way, if the MATLAB routine FGETS is used, instead of
%    FGETL, then the variable LINE will include line termination
%    characters, which means that a blank line would not actually
%    have zero characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Output, integer ROW_NUM, the number of rows found.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_ROW_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_ROW_COUNT - Error!' );
  end

  blank_num = 0;
  comment_num = 0;
  row_num = 0;

  record_num = 0;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    record_num = record_num + 1;
    record_length = s_len_trim ( line );

    if ( record_length <= 0 )
      blank_num = blank_num + 1;
    elseif ( line(1) == '#' )
      comment_num = comment_num + 1;
    else
      row_num = row_num + 1;
    end

  end

  fclose ( input_unit );

  return
end
function isgn = i4col_compare ( m, n, a, i, j )

%*****************************************************************************80
%
%% I4COL_COMPARE compares columns I and J of a integer array.
%
%  Example:
%
%    Input:
%
%      M = 3, N = 4, I = 2, J = 4
%
%      A = (
%        1  2  3  4
%        5  6  7  8
%        9 10 11 12 )
%
%    Output:
%
%      ISGN = -1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, integer A(M,N), an array of N columns of vectors of length M.
%
%    Input, integer I, J, the columns to be compared.
%    I and J must be between 1 and N.
%
%    Output, integer ISGN, the results of the comparison:
%    -1, column I < column J,
%     0, column I = column J,
%    +1, column J < column I.
%

%
%  Check.
%
  if ( i < 1)
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4COL_COMPARE - Fatal error!\n' );
    fprintf ( 1, '  Column index I = %d < 1.\n', i );
    error ( 'I4COL_COMPARE - Fatal error!' );
  end

  if ( n < i )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4COL_COMPARE - Fatal error!\n' );
    fprintf ( 1, '  N = %d < column index I = %d.\n', n, i );
    error ( 'I4COL_COMPARE - Fatal error!' );
  end

  if ( j < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4COL_COMPARE - Fatal error!\n' );
    fprintf ( 1, '  Column index J = %d < 1.\n', j );
    error ( 'I4COL_COMPARE - Fatal error!' );
  end

  if ( n < j )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4COL_COMPARE - Fatal error!\n' );
    fprintf ( 1, '  N = %d < column index J = %d.\n', n, j );
    error ( 'I4COL_COMPARE - Fatal error!' );
  end

  isgn = 0;

  if ( i == j )
    return
  end

  k = 1;

  while ( k <= m )

    if ( a(k,i) < a(k,j) )
      isgn = -1;
      return
    elseif ( a(k,j) < a(k,i) )
      isgn = +1;
      return
    end

    k = k + 1;

  end

  return
end
function a = i4col_sort_a ( m, n, a )

%*****************************************************************************80
%
%% I4COL_SORT_A ascending sorts an I4COL.
%
%  Discussion:
%
%    In lexicographic order, the statement "X < Y", applied to two real
%    vectors X and Y of length M, means that there is some index I, with
%    1 <= I <= M, with the property that
%
%      X(J) = Y(J) for J < I,
%    and
%      X(I) < Y(I).
%
%    In other words, the first time they differ, X is smaller.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows of A, and the length of
%    a vector of data.
%
%    Input, integer N, the number of columns of A.
%
%    Input, integer A(M,N), the array of N columns of M-vectors.
%
%    Output, integer A(M,N), the columns of A have been sorted in ascending
%    lexicographic order.
%
  if ( m <= 0 )
    return
  end

  if ( n <= 1 )
    return
  end
%
%  Initialize.
%
  indx = 0;
  isgn = 0;
%
%  Call the external heap sorter.
%
  while ( 1 )

    [ indx, i, j ] = sort_heap_external ( n, indx, isgn );
%
%  Interchange the I and J objects.
%
    if ( 0 < indx )

      a = i4col_swap ( m, n, a, i, j );
%
%  Compare the I and J objects.
%
    elseif ( indx < 0 )

      isgn = i4col_compare ( m, n, a, i, j );

    elseif ( indx == 0 )

      break

    end

  end

  return
end
function a = i4col_swap ( m, n, a, i, j )

%*****************************************************************************80
%
%% I4COL_SWAP swaps columns I and J of a integer array of column data.
%
%  Example:
%
%    Input:
%
%      M = 3, N = 4, I = 2, J = 4
%
%      A = (
%        1  2  3  4
%        5  6  7  8
%        9 10 11 12 )
%
%    Output:
%
%      A = (
%        1  4  3  2
%        5  8  7  6
%        9 12 11 10 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns in the array.
%
%    Input, integer A(M,N), an array of N columns of length M.
%
%    Input, integer I, J, the columns to be swapped.
%
%    Output, integer A(M,N), the array, with columns I and J swapped.
%
  if ( i < 1 || n < i || j < 1 || n < j )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4COL_SWAP - Fatal error!\n' );
    fprintf ( 1, '  I or J is out of bounds.\n' );
    fprintf ( 1, '  I =    %d\n', i );
    fprintf ( 1, '  J =    %d\n', j );
    fprintf ( 1, '  N =    %d\n', n );
    error ( 'I4COL_SWAP - Fatal error!' );
  end

  if ( i == j )
    return
  end

  col(1:m) = a(1:m,i)';
  a(1:m,i) = a(1:m,j);
  a(1:m,j) = col(1:m)';

  return
end
function [ j1, j2, j3 ] = i4i4i4_sort_a ( i1, i2, i3 )

%*****************************************************************************80
%
%% I4I4I4_SORT_A ascending sorts a triple of integers.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I1, I2, I3, the values to sort.
%
%    Output, integer J1, J2, J3, the sorted values.
%
  j1 = min ( min ( i1, i2 ), min ( i2, i3 ) );
  j2 = min ( max ( i1, i2 ), ...
       min ( max ( i2, i3 ), max ( i3, i1 ) ) );
  j3 = max ( max ( i1, i2 ), max ( i2, i3 ) );

  return
end
function table = i4mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% I4MAT_DATA_READ reads data from an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns in the data.
%
%    Output, integer TABLE(M,N), the point coordinates.
%
  table = zeros ( m, n );
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %d' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the input file.\n' );
    error ( 'I4MAT_DATA_READ - Error!' );
  end

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'I4MAT_DATA_READ - Error!\n' );
      fprintf ( 1, '  End of input while reading data.\n' );
      error ( 'I4MAT_DATA_READ - Error!' );
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = i4mat_header_read ( input_filename )

%*****************************************************************************80
%
%% I4MAT_HEADER_READ reads the header from an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

  return
end
function i4mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% I4MAT_TRANSPOSE_PRINT_SOME prints some of an I4MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, integer A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, an optional title.
%
  incx = 10;

  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d  ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );
    fprintf ( 1, '\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d  ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%7d  ', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end
function i4mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% I4MAT_WRITE writes an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 July 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer TABLE(M,N), the points.
%
%    Input, logical HEADER, is TRUE if the header is to be included.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'I4MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %12d', round ( table(i,j) ) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function element_node = mesh_base_one ( element_order, element_num, ...
  element_node )

%*****************************************************************************80
%
%% MESH_BASE_ONE ensures that the element definition is one-based.
%
%  Discussion:
%
%    The ELEMENT_NODE array contains nodes indices that form elements.
%    The convention for node indexing might start at 0 or at 1.
%    Since a MATLAB program will naturally assume a 1-based indexing, it is
%    necessary to check a given element definition and, if it is actually
%    0-based, to convert it.
%
%    This function attempts to detect 0-based node indexing and correct it.
%
%    Thanks to Feifei Xu for pointing out that I was subtracting 1 when I
%    should have been adding 1!  29 November 2012.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 November 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ELEMENT_ORDER, the order of the elements.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input/output, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM), the element
%    definitions.
%
  node_min = min ( min ( element_node(1:element_order,1:element_num) ) );

  if ( node_min == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MESH_BASE_ONE:\n' );
    fprintf ( 1, '  The element indexing appears to be 0-based!\n' );
    fprintf ( 1, '  This will be converted to 1-based.\n' );
    element_node(1:element_order,1:element_num) = ...
      element_node(1:element_order,1:element_num) + 1;
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MESH_BASE_ONE:\n' );
    fprintf ( 1, '  The element indexing appears to be 1-based!\n' );
    fprintf ( 1, '  No conversion is necessary.\n' );
  end

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
%% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LEN, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end
function word_num = s_word_count ( s )

%*****************************************************************************80
%
%% S_WORD_COUNT counts the number of "words" in a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be examined.
%
%    Output, integer WORD_NUM, the number of "words" in the string.
%    Words are presumed to be separated by one or more blanks.
%
  FALSE = 0;
  TRUE = 1;

  word_num = 0;
  s_length = length ( s );

  if ( s_length <= 0 )
    return;
  end

  blank = TRUE;

  for i = 1 : s_length

    if ( s(i) == ' ' )
      blank = TRUE;
    elseif ( blank == TRUE )
      word_num = word_num + 1;
      blank = FALSE;
    end

  end

  return
end
function [ indx, i, j ] = sort_heap_external ( n, indx, isgn )

%*****************************************************************************80
%
%% SORT_HEAP_EXTERNAL externally sorts a list of items into ascending order.
%
%  Discussion:
%
%    The actual list of data is not passed to the routine.  Hence this
%    routine may be used to sort integers, reals, numbers, names,
%    dates, shoe sizes, and so on.  After each call, the routine asks
%    the user to compare or interchange two items, until a special
%    return value signals that the sorting is completed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 February 2004
%
%  Author:
%
%    Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf.
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Albert Nijenhuis, Herbert Wilf.
%    Combinatorial Algorithms,
%    Academic Press, 1978, second edition,
%    ISBN 0-12-519260-6.
%
%  Parameters:
%
%    Input, integer N, the number of items to be sorted.
%
%    Input, integer INDX, the main communication signal.
%    The user must set INDX to 0 before the first call.
%    Thereafter, the user should set the input value of INDX
%    to the output value from the previous call.
%
%    Input, integer ISGN, results of comparison of elements I and J.
%    (Used only when the previous call returned INDX less than 0).
%    ISGN <= 0 means I is less than or equal to J;
%    0 <= ISGN means I is greater than or equal to J.
%
%    Output, integer INDX, the main communication signal.
%    If INDX is
%
%      greater than 0, the user should:
%      * interchange items I and J;
%      * call again.
%
%      less than 0, the user should:
%      * compare items I and J;
%      * set ISGN = -1 if I < J, ISGN = +1 if J < I;
%      * call again.
%
%      equal to 0, the sorting is done.
%
%    Output, integer I, J, the indices of two items.
%    On return with INDX positive, elements I and J should be interchanged.
%    On return with INDX negative, elements I and J should be compared, and
%    the result reported in ISGN on the next call.
%
  persistent i_save;
  persistent j_save;
  persistent k;
  persistent k1;
  persistent n1;
%
%  INDX = 0: This is the first call.
%
  if ( indx == 0 )

    i_save = -1;
    j_save = -1;
    k = floor ( n / 2 );
    k1 = k;
    n1 = n;
%
%  INDX < 0: The user is returning the results of a comparison.
%
  elseif ( indx < 0 )

    if ( indx == -2 )

      if ( isgn < 0 )
        i_save = i_save + 1;
      end

      j_save = k1;
      k1 = i_save;
      indx = -1;
      i = i_save;
      j = j_save;
      return;
    end

    if ( 0 < isgn )
      indx = 2;
      i = i_save;
      j = j_save;
      return;
    end

    if ( k <= 1 )

      if ( n1 == 1 )
        i_save = 0;
        j_save = 0;
        indx = 0;
      else
        i_save = n1;
        n1 = n1 - 1;
        j_save = 1;
        indx = 1;
      end

      i = i_save;
      j = j_save;
      return;

    end

    k = k - 1;
    k1 = k;
%
%  0 < INDX, the user was asked to make an interchange.
%
  elseif ( indx == 1 )

    k1 = k;

  end

  while ( 1 )

    i_save = 2 * k1;

    if ( i_save == n1 )
      j_save = k1;
      k1 = i_save;
      indx = -1;
      i = i_save;
      j = j_save;
      return;
    elseif ( i_save <= n1 )
      j_save = i_save + 1;
      indx = -2;
      i = i_save;
      j = j_save;
      return;
    end

    if ( k <= 1 )
      break;
    end

    k = k - 1;
    k1 = k;

  end

  if ( n1 == 1 )
    i_save = 0;
    j_save = 0;
    indx = 0;
    i = i_save;
    j = j_save;
  else
    i_save = n1;
    n1 = n1 - 1;
    j_save = 1;
    indx = 1;
    i = i_save;
    j = j_save;
  end

  return
end
function tetra_neighbor = tet_mesh_neighbor_tets ( tetra_order, tetra_num, ...
  tetra_node )

%*****************************************************************************80
%
% TET_MESH_NEIGHBOR_TETS determines tetrahedron neighbors.
%
%  Discussion:
%
%    A tet mesh of a set of nodes can be completely described by
%    the coordinates of the nodes, and the list of nodes that make up
%    each tetrahedron.  In the most common case, four nodes are used.
%    There is also a 10 node case, where nodes are also placed on
%    the midsides of the tetrahedral edges.
%
%    This routine can handle 4 or 10-node tetrahedral meshes.  The
%    10-node case is handled simply by ignoring the six midside nodes,
%    which are presumed to be listed after the vertices.
%
%    The tetrahedron adjacency information records which tetrahedron
%    is adjacent to a given tetrahedron on a particular face.
%
%    This routine creates a data structure recording this information.
%
%    The primary amount of work occurs in sorting a list of 4 * TETRA_NUM
%    data items.
%
%    The neighbor tetrahedrons are indexed by the face they share with
%    the tetrahedron.
%
%    Each face of the tetrahedron is indexed by the node which is NOT
%    part of the face.  That is:
%
%    * Neighbor 1 shares face 1 defined by nodes 2, 3, 4.
%    * Neighbor 2 shares face 2 defined by nodes 1, 3, 4;
%    * Neighbor 3 shares face 3 defined by nodes 1, 2, 4;
%    * Neighbor 4 shares face 4 defined by nodes 1, 2, 3.
%
%    For instance, if the (transposed) TETRA_NODE array was:
%
%    Row       1      2      3      4
%    Col
%
%      1       4      3      5      1
%      2       4      2      5      1
%      3       4      7      3      5
%      4       4      7      8      5
%      5       4      6      2      5
%      6       4      6      8      5
%
%    then the (transposed) TETRA_NEIGHBOR array should be:
%
%    Row       1      2      3      4
%    Col
%
%      1      -1      2     -1      3
%      2      -1      1     -1      5
%      3      -1      1      4     -1
%      4      -1      6      3     -1
%      5      -1      2      6     -1
%      6      -1      4      5     -1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer TETRA_ORDER, the order of the tetrahedrons.
%
%    Input, integer TETRA_NUM, the number of tetrahedrons.
%
%    Input, integer TETRA_NODE(TETRA_ORDER,TETRA_NUM), the nodes that make up
%    each tetrahedron.
%
%    Output, integer TETRA_NEIGHBOR(4,TETRA_NUM), the four tetrahedrons that
%    are direct neighbors of a given tetrahedron.  If there is no neighbor
%    sharing a given face, the index is set to -1.
%

%
%  Step 1.
%  From the list of nodes for tetrahedron T, of the form: (I,J,K,L)
%  construct the four face relations:
%
%    (J,K,L,1,T)
%    (I,K,L,2,T)
%    (I,J,L,3,T)
%    (I,J,K,4,T)
%
%  In order to make matching easier, we reorder each triple of nodes
%  into ascending order.
%
  faces = zeros ( 5, 4 * tetra_num );
  for tetra = 1 : tetra_num

    i = tetra_node(1,tetra);
    j = tetra_node(2,tetra);
    k = tetra_node(3,tetra);
    l = tetra_node(4,tetra);

    [ a, b, c ] = i4i4i4_sort_a ( j, k, l );

    faces(1:5,4*(tetra-1)+1) = [ a, b, c, 1, tetra ]';

    [ a, b, c ] = i4i4i4_sort_a ( i, k, l );

    faces(1:5,4*(tetra-1)+2) = [ a, b, c, 2, tetra ]';

    [ a, b, c ] = i4i4i4_sort_a ( i, j, l );

    faces(1:5,4*(tetra-1)+3) = [ a, b, c, 3, tetra ]';

    [ a, b, c ] = i4i4i4_sort_a ( i, j, k );

    faces(1:5,4*(tetra-1)+4) = [ a, b, c, 4, tetra ]';

  end
%
%  Step 2. Perform an ascending dictionary sort on the neighbor relations.
%  We only intend to sort on rows 1:3; the routine we call here
%  sorts on rows 1 through 5 but that won't hurt us.
%
%  What we need is to find cases where two tetrahedrons share a face.
%  By sorting the columns of the FACES array, we will put shared faces
%  next to each other.
%
  faces = i4col_sort_a ( 5, 4*tetra_num, faces );
%
%  Step 3. Neighboring tetrahedrons show up as consecutive columns with
%  identical first three entries.  Whenever you spot this happening,
%  make the appropriate entries in TETRA_NEIGHBOR.
%
  tetra_neighbor(1:4,1:tetra_num) = -1;

  face = 1;

  while ( 1 )

    if ( 4 * tetra_num <= face )
      break
    end

    if ( faces(1:3,face) == faces(1:3,face+1) )
      face1 = faces(4,face);
      tetra1 = faces(5,face);
      face2 = faces(4,face+1);
      tetra2 = faces(5,face+1);
      tetra_neighbor(face1,tetra1) = tetra2;
      tetra_neighbor(face2,tetra2) = tetra1;
      face = face + 2;
    else
      face = face + 1;
    end

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
clc;
close all;
clearvars;
rng('default');
n = 10;
% Let's get a symmetric positive definite matrix
A = gallery('gcdmat', n);
% Compute its Cholesky decomposition
L = chol(A, 'lower')
b = randn(n, 1);
% solve the problem
x = L \b
x = spx.la.tris.forward_col(L, b)

author: indigits

repo: sparse-plex


text:
function [mx, mt] = vl_nnmasknan(x, t, varargin)
%VL_NNMASKNAN masks out NaN data from targets
%   [MX, MT] = VL_NNMASKNAN(X,T) masks out NaN data present in targets
%   T and drops the corresponding entries in X.  Produces masked outputs
%   MX and MT where MT consists of the non-NaN entries of T and MX consists
%   of the corresponding entries of X. The shape of both MX and MT will be
%   Nx1, where N is the number of non-NaN entries of T.
%
%   NOTE: I think this function is borrowed from somewhere, since I don't
%   remember writing it.  Unfortunately I can't remember where it might
%   be from, making attribution difficult.
%
% Copyright (C) 2017 Samuel Albanie
% Licensed under The MIT License [see LICENSE.md for details]

  [~, dzdy] = vl_argparsepos(struct(), varargin) ;

  mask = isnan(t) ;
  mx = x ; mt = t ;

  if isempty(dzdy)
    mx(mask(:)) = [] ;
    mt(mask(:)) = [] ;
    mx = mx' ; mt = mt' ;
  else
    mx = zeros(size(x), 'like', x) ;
    mx(~mask) = dzdy{1} ;
  end

repo: MFT


text:
%% simple example of use of renorm_sibling_3d_same_scale
clear; close all;
x = uiuc_sample;
%% compute roto-translation scattering of an image
options.Q = 1;
options.J = 5;
Wop = wavelet_factory_3d_pyramid(options, options, options);
Sx = scat(x, Wop);

%% L1 renormalization
op = @(x)(sum(x, 3));
Sx_renorm = renorm_sibling_3d_same_scale(Sx, op);

%% L1 + smoothing renormalization
op = renorm_factory_L1_smoothing(2);
[Sx_renorm, siblings] = renorm_sibling_3d_same_scale(Sx, op);

%%
image_scat(Sx, 0, 0);
%%
image_scat(Sx_renorm, 0, 0);
%%
close all;
imagesc(min(image_scat_layer(Sx_renorm{3},0,0),0.2));

%%
image_scat(Sx_renorm, 1, 1);
%%
image_scat(Sx, 1, 1);

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
% test for farthest point sampling for volumetric mesh

n = 150;

name = 'constant';
name = 'central';
name = 'sphere';
name = 'medical';

rep = ['results/farthest-sampling-3d/' name '/'];
if exist(rep)~=7
    mkdir(rep);
end

switch name
    case 'constant'
        W = ones(n,n,n);
    case {'central' 'sphere'}
        x = linspace(-1,1,n);
        [X,Y,Z] = meshgrid(x, x, x);
        W = sqrt( X.^2+Y.^2+Z.^2 );
        if strcmp(name, 'sphere')
            r = .6;
            W = rescale( -(W<=r) );
        end
        W = rescale(W,.1,1);
    case 'medical'
        load brain1-crop-256.mat
        M = crop(M,n);
        W = rescale(-abs(M-median(M(:))));
        W = rescale( clamp(W,.6,1),.2,1);
    otherwise
        error('Unknown potential');
end

method = 'geodesic';
method = 'delaunay';

% padd to treat correctly boundary
if strcmp(method, 'geodesic')
    k = 10;
    W1 = perform_image_extension(W,n+2*k, '2side');
end

% plot sampling location
ms = 20; lw = 3;
i = 0;
for npoints = [400 800 1600 3000 5000 10000]
    i = i+1;
    disp('Perform farthest point sampling');
    if i==1
        [x,y,z] = meshgrid([1 n], [1 n], [1 n]);
        vertex = cat(1, x(:)', y(:)', z(:)');
    end
    vertex = perform_farthest_point_sampling( W, vertex, npoints-size(vertex,2) );

    % display sampling
    clf;
    hold on;
    h = imageplot(W);
    view(60,30); zoom(.8);  camlight;
    alphamap(linspace(0,.1,100)); vol3d(h);
    h = plot3(vertex(1,:),vertex(2,:),vertex(3,:), '.');
    set(h, 'MarkerSize', ms);
    hold off;
    saveas(gcf, [rep name '-sampling-' num2string_fixeddigit(npoints,5) '.png'], 'png');

    % compute the associated triangulation
    if strcmp(method, 'delaunay')
        faces = delaunay3(vertex(1,:),vertex(2,:),vertex(3,:))';
    else
        vertex1 = vertex+k;
        [D,Z,Q] = perform_fast_marching(W1, vertex1);
        faces = compute_voronoi_triangulation(Q,vertex);
    end

    w = [1 1 1]; w = w(:)/sqrt(sum(w.^2));
    t = sum(vertex.*repmat(w,[1 size(vertex,2)]));
    delta = max(t)-min(t);
    offlist = linspace( min(t(:)) + .3*delta, max(t(:)) - .3*delta, 5 );
    for i=1:length(offlist)
        options.cutting_plane = w;
        options.cutting_offs = offlist(i);
        clf;
        hold on;
        %h = imageplot(W);
        %alphamap(linspace(0,.1,100));
        plot_mesh(vertex,faces, options);
        % set(h{2}, 'Marker', 'r.');
        view(60,30); zoom(.8);  camlight;
        % vol3d(h);
        hold off;
        saveas(gcf, [rep name '-sampling-' num2string_fixeddigit(npoints,5) '-mesh-' num2str(i) '.png'], 'png');
    end

end

author: alecjacobson

repo: gptoolbox


text:
function varargout = mrmMakeMovieGUI(varargin)
% mrmMakeMovieGUI(meshID)
%   Given a mesh ID #, opens a GUI allowing you to set up a series of
%   events to be displayed in a .avi file.
%
%   Events consist of waypoints, transitions, and pauses.  All pauses must
%   follow waypoints.  Transitions must appear with waypoints on either
%   side.  Move the mesh manually to find the desired waypoint location, or
%   use the presets.
%
% Usage:
%   mrmMakeMovieGUI(meshID);
%
% Known issues:
%   Jump to view isn't accurate with diffusion meshes for whatever reason.
%   Probably need to edit mrmRotateCamera and add a conditional for when
%   we're rotating such meshes, providing different coordinates (should
%   they be consistent across diffusion meshes).
%
% [renobowen@gmail.com 2010]
%

    % Begin initialization code - DO NOT EDIT
    gui_Singleton = 1;
    gui_State = struct('gui_Name',       mfilename, ...
                       'gui_Singleton',  gui_Singleton, ...
                       'gui_OpeningFcn', @mrmMakeMovieGUI_OpeningFcn, ...
                       'gui_OutputFcn',  @mrmMakeMovieGUI_OutputFcn, ...
                       'gui_LayoutFcn',  [] , ...
                       'gui_Callback',   []);
    if nargin && ischar(varargin{1})
        gui_State.gui_Callback = str2func(varargin{1});
    end

    if nargout
        [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
    else
        gui_mainfcn(gui_State, varargin{:});
    end
    % End initialization code - DO NOT EDIT

end

function mrmMakeMovieGUI_OpeningFcn(hObject, eventdata, handles, meshID)
    if (~exist('meshID', 'var')), meshID = []; end

    % Choose default command line output for mrmMakeMovieGUI
    handles.output      = hObject;
    handles.meshID      = meshID;
    handles.events      = {};
    handles.isSaved     = true;
    handles.filename    = [];
    handles.filedir     = [];
    handles.clipboard   = [];
    handles.index       = 0;

    set(handles.EventTypeMenu, 'String', {'Waypoint', 'Transition', 'Pause'});

    % Update handles structure
    handles = UpdateFileStatus(handles, true);
    guidata(hObject, handles);

end

function varargout = mrmMakeMovieGUI_OutputFcn(hObject, eventdata, handles)
    if (isempty(handles.meshID))
        fprintf(1, 'No mesh ID # specified - closing GUI.\n');
        close(handles.MakeMovieGUI);
    end
end

function handles = TurnOffInvalidButtons(handles)
    if (isempty(handles.clipboard))
        enablePaste = 'off';
    else
        enablePaste = 'on';
    end
    if (length(handles.events) < 1)
        enableEventDependent = 'off';
    else
        enableEventDependent = 'on';
    end
    % Paste buttons
    set(handles.EventMenu_Paste, 'Enable', enablePaste);
    set(handles.EventContextMenu_Paste, 'Enable', enablePaste);

    % Buttons dependent on existence of events
    set(handles.ExportButton, 'Enable', enableEventDependent);
    set(handles.EditEventButton, 'Enable', enableEventDependent);

    set(handles.PreviewButton, 'Enable', enableEventDependent);
    set(handles.DeleteEventButton, 'Enable', enableEventDependent);

    set(handles.MoveDownButton, 'Enable', enableEventDependent);
    set(handles.MoveUpButton, 'Enable', enableEventDependent);

    set(handles.MovieMenu_Export, 'Enable', enableEventDependent);
    set(handles.MovieMenu_Preview, 'Enable', enableEventDependent);
    set(handles.EventMenu_Delete, 'Enable', enableEventDependent);
    set(handles.EventMenu_Edit, 'Enable', enableEventDependent);

    set(handles.EventMenu_Copy, 'Enable', enableEventDependent);
    set(handles.EventMenu_Cut, 'Enable', enableEventDependent);

    set(handles.EventContextMenu_Delete, 'Enable', enableEventDependent);
    set(handles.EventContextMenu_Edit, 'Enable', enableEventDependent);

    set(handles.EventContextMenu_Copy, 'Enable', enableEventDependent);
    set(handles.EventContextMenu_Cut, 'Enable', enableEventDependent);

    set(handles.FileMenu_SaveAs, 'Enable', enableEventDependent);
    set(handles.FileMenu_Save, 'Enable', enableEventDependent);
end

function handles = UpdateFileStatus(handles, saving)
    if (~isempty(handles.filename))
        if (saving)
            handles.isSaved = true;
            set(handles.FilenameText, 'String', handles.filename);
        else
            handles.isSaved = false;
            set(handles.FilenameText, 'String', [handles.filename '*']);
        end
    else
        if (saving)
            handles.isSaved = true;
        else
            handles.isSaved = false;
        end
        set(handles.FilenameText, 'String', 'No save file.');
    end
    handles = TurnOffInvalidButtons(handles);
end

function EventList_Callback(hObject, eventdata, handles)
    if (length(handles.events) < 1), return; end

    index = get(handles.EventList, 'Value');
    event = handles.events{index};
    if (strcmp(event.eventType,'waypoint'))
        mrmRotateCamera(handles.meshID, event.rotation, [], event.frustum, [], event.origin);
    end
end

function MoveUpButton_Callback(hObject, eventdata, handles)
    index = get(handles.EventList, 'Value');

    ExchangeEvents(hObject, handles, index, index - 1);

end

function MoveDownButton_Callback(hObject, eventdata, handles)
    index = get(handles.EventList, 'Value');

    ExchangeEvents(hObject, handles, index, index + 1);

end

function ExchangeEvents(hObject, handles, selected, exchanged)
    strings = get(handles.EventList, 'String');
    if ((exchanged > length(strings)) || exchanged < 1), return; end

    tmpString           = strings{exchanged};
    strings{exchanged}  = strings{selected};
    strings{selected}   = tmpString;

    tmpEvent                    = handles.events{exchanged};
    handles.events{exchanged}   = handles.events{selected};
    handles.events{selected}    = tmpEvent;

    selected = exchanged;
    set(handles.EventList, 'String', strings);
    set(handles.EventList, 'Value', selected);

	handles = UpdateFileStatus(handles, false);
    guidata(hObject, handles);

end

function AddEvent_Callback(hObject, eventdata, handles, eventType, doInsert)
    % Navigating around a godawful bug with single selection GUIs
    if (length(handles.events) < 1), set(handles.EventList, 'Value', 1); end

    % Get event type from drop down if it's not specified
    if (~exist('eventType', 'var') || isempty(eventType))
        eventTypes = get(handles.EventTypeMenu, 'String');
        eventType = eventTypes{get(handles.EventTypeMenu, 'Value')};
    end

    % Ship off task of getting parameters about event to smaller GUIs
    switch lower(eventType)
        case 'waypoint'
            waitfor(mrmMakeMovieGUI_Waypoint(handles.MakeMovieGUI, handles.meshID));

        case 'transition'
            waitfor(mrmMakeMovieGUI_Transition(handles.MakeMovieGUI));

        case 'pause'
            waitfor(mrmMakeMovieGUI_Pause(handles.MakeMovieGUI));

        otherwise
            fprintf(1, 'Unrecognized radio button selected.\n');

    end

    % Get data back from the GUIs, if any
    data = get(handles.MakeMovieGUI, 'UserData'); % Get the data back from the pop up
    if (isempty(data)), return; end % We're done if there are no params

    % If we're not inserting, tack it after current selection
    numEvents = length(handles.events);
    if (~exist('doInsert', 'var') || isempty(doInsert) || (numEvents < 1))
        index = numEvents + 1;
    else
        index = get(handles.EventList, 'Value') + 1;
    end

    handles = InsertEvent(hObject, handles, data, index);

    set(handles.MakeMovieGUI, 'UserData', []); % Done with the GUI data, clear it
    guidata(hObject, handles);
end

function handles = InsertEvent(hObject, handles, event, position)
% Insert event into visible list and storage array at specified position
    handles.events = CellInsert(handles.events, event, position);
    strings = get(handles.EventList, 'String');
    strings = CellInsert(strings, GetListString(event), position);
    set(handles.EventList, 'Value', position);
    set(handles.EventList, 'String', strings);
    handles = UpdateFileStatus(handles, false);
end

function cells = CellInsert(cells, value, index)
% Insert a value into a cell, shifting over the values in its place if
% necessary
    numCells = length(cells);
    if (numCells >= index)
        for i = numCells:-1:index
            cells{i+1} = cells{i};
        end
    end
    cells{index} = value;
end

function listString = GetListString(data)
    switch (data.eventType)
        case 'waypoint'
            string = '[WAYPOINT]';

        case 'transition'
            string = '[TRANSITION]';

        case 'pause'
            string = '[PAUSE]';

        otherwise
            fprintf(1, 'Unrecognized event type ''%s''', data.eventType);
            listString = [];
            return;

    end

    listString = sprintf('%s %s', string, data.label);
end

function DeleteEvent_Callback(hObject, eventdata, handles)
    index = get(handles.EventList, 'Value');

    strings = get(handles.EventList, 'String');

    strings{index}              = [];
    handles.events{index}       = [];

    strings             = cellRemoveEmpty(strings);
    handles.events      = cellRemoveEmpty(handles.events);

    set(handles.EventList, 'String', strings);

    numEvents = length(handles.events);
    if (numEvents == 0)
        index = 1;
    else
        if (index > numEvents)
            index = numEvents;
        end
    end

    set(handles.EventList, 'Value', index);

	handles = UpdateFileStatus(handles, false);
    guidata(hObject, handles);

end

function frames = CountFrames(events)
    frames = 0;
    numEvents = length(events);
    for i = 1:numEvents
        event = events{i};
        if (isfield(event, 'frames'))
            frames = frames + event.frames;
        end
    end
end

function RenderMovie_Callback(hObject, eventdata, handles, preview)
    if (AreEventsWellFormed(handles))
        frames = CountFrames(handles.events);
        framesPerSec    = str2num(get(handles.FPSTextField, 'String'));
        events          = handles.events;
        nEvents         = length(events);
        lastWaypoint    = [];

        if (~preview) % Set up for saving out movie
            [movFile, movDir]   = uiputfile('brainMovie.avi');
            if (movFile == 0), return; end

        	host                = 'localhost';
            f.filename          = 'nosave';
            frame               = 1;
            ShowProgress(handles, movFile, 0);
            % Skipping preallocation of M, will do so if it's too slow
        end

        for i = 1:nEvents
            event = handles.events{i};
            switch (event.eventType)
                case 'transition'
                    [rotation frustum origin] = BuildTransition(event, lastWaypoint, GetNextWaypoint(events, i));
                    for j = 1:event.frames
                        mrmRotateCamera(handles.meshID, rotation(j, :), [], frustum(j, :), [], origin(j, :));
                        if (~preview)
                            [id,stat,res] = mrMesh(host, handles.meshID, 'screenshot', f);
                            ShowProgress(handles, movFile, frame/frames * 100);
                            M(frame) = im2frame(permute(res.rgb, [2,1,3])./255);
                            frame = frame + 1;
                        end
                    end
                case 'waypoint'
                    lastWaypoint = event;
                    mrmRotateCamera(handles.meshID, event.rotation, [], event.frustum, [], event.origin);
                case 'pause'
                    if (preview)
                        pause(event.frames/framesPerSec);
                    else
                        [id,stat,res] = mrMesh(host, handles.meshID, 'screenshot', f);
                        for j = 1:event.frames
                            ShowProgress(handles, movFile, frame/frames * 100);
                            M(frame) = im2frame(permute(res.rgb, [2,1,3])./255);
                            frame = frame + 1;
                        end
                    end
                otherwise
                    fprintf(1, 'Unrecognized event type ''%s''', data.eventType);
                    return;
            end
        end

        if (~preview)
            if ispc
                % Compressors for Windows are annoying.  None works.  Not
                % sure if the others in doc movie2avi work.
                movie2avi(M, fullfile(movDir,movFile), ...
                    'FPS', framesPerSec, ...
                    'Compression','RLE');
            else
                % Unix/Mac
                movie2avi(M, fullfile(movDir,movFile), ...
                    'FPS', framesPerSec, ...
                    'Compression','None');

            end

            set(handles.ProgressText, 'String', []);
        end

    else
        errordlg('Invalid events list.  See wiki for guidelines.');
    end

end

function ShowProgress(handles, filename, percent)
    set(handles.ProgressText, 'String', sprintf('Exporting %s... %2.0f%%', filename, percent));
end

function [rotation frustum origin] = BuildTransition(transition, waypointStart, waypointEnd)
    nFrames     = transition.frames;
    rotation    = zeros(nFrames, 3);
    rotStart    = waypointStart.rotation;
    rotEnd      = waypointEnd.rotation;

    for i = 1:3
        if (rotEnd(i) == rotStart(i))
            rotation(:, i) = rotStart(i);
            continue;
        end

        switch (transition.rotate{i})
            case '+'
                if (rotEnd(i) < rotStart(i))
                    rotEnd(i) = pi + (pi + rotEnd(i));
                end

            case '-'
                if (rotEnd(i) > rotStart(i))
                    rotEnd(i) = -pi - (pi - rotEnd(i));
                end

            case 'FIX'
                rotation(:,i) = rotStart(i);
                continue;

            otherwise
                fprintf(1, 'Unrecognized rotation direction ''%s''', transition.rotate{i});
        end

        if (rotEnd(i) == rotStart(i))
            rotation(:, i) = rotStart(i);
            continue;
        end

        delta           = rotEnd(i) - rotStart(i);
        rotation(:, i)  = rotStart(i):(delta/(nFrames - 1)):rotEnd(i);
    end

    frustum = InterpolateFrames(waypointStart.frustum, waypointEnd.frustum, nFrames);
    origin  = InterpolateFrames(waypointStart.origin, waypointEnd.origin, nFrames);

end

function vector = InterpolateFrames(start, finish, nFrames)
    nEntries = length(start);
    vector = zeros(nFrames, nEntries);
    for i = 1:nEntries
        if (start(i) == finish(i))
            vector(:,i) = start(i);
            continue;
        end

        delta       = finish(i) - start(i);
        vector(:,i) = start(i):(delta/(nFrames - 1)):finish(i);
    end

end

function waypoint = GetNextWaypoint(events, startInd)
    waypoint = [];
    for i = startInd:length(events)
       if (strcmp(events{i}.eventType, 'waypoint'))
           waypoint = events{i};
           return;
       end
    end
end

function bool = AreEventsWellFormed(handles)
    events  = handles.events;
    nEvents = length(events);

    if (nEvents < 2) % Need at least 2 events (waypoint + pause) to proceed
        bool = false;
        return;
    end

    event = events{1};
    if (~(strcmp(event.eventType, 'waypoint'))) % First event needs to be a waypoint
        bool = false;
        return;
    end

    balanceCounter  = 0;
    hasFrames       = 0;
    for i = 2:nEvents
        event = events{i};
        if (~hasFrames && (strcmp(event.eventType, 'pause') || strcmp(event.eventType, 'transition')))
            hasFrames = 1;
        end

        if (strcmp(event.eventType, 'transition'))
            if (balanceCounter), bool = false; return; end
            balanceCounter = balanceCounter + 1;
        elseif (strcmp(event.eventType, 'waypoint'))
            if (~balanceCounter), continue; end
            balanceCounter = balanceCounter - 1;
        elseif (strcmp(event.eventType, 'pause'))
            if (balanceCounter), bool = false; return; end
        end
    end
    if (~hasFrames || balanceCounter), bool = false; return; end
    bool = true;
end

function EditEvent_Callback(hObject, eventdata, handles)
    events = handles.events;
    nEvents = length(events);
    if (nEvents < 1), return; end
    selected = get(handles.EventList, 'Value');
    editData = events{selected};

    switch (editData.eventType)
        case 'waypoint'
            waitfor(mrmMakeMovieGUI_Waypoint(handles.MakeMovieGUI, handles.meshID, editData));

        case 'transition'
            waitfor(mrmMakeMovieGUI_Transition(handles.MakeMovieGUI, editData));

        case 'pause'
            waitfor(mrmMakeMovieGUI_Pause(handles.MakeMovieGUI, editData));

        otherwise
            fprintf(1, 'Unrecognized radio button selected.\n');

    end

    editData = get(handles.MakeMovieGUI, 'UserData');
    if (isempty(editData)), return; end

    handles.events{selected} = editData;

    strings = get(handles.EventList, 'String');
    strings{selected} = GetListString(editData);
    set(handles.EventList, 'String', strings);

    handles = UpdateFileStatus(handles, false);
    set(handles.MakeMovieGUI, 'UserData', []); % Clear it out to not muck up future use of this var
    guidata(hObject, handles);
end

function SaveData(handles)
    if (~isempty(handles.filename))
        fileID      = 'EventsList';
        version     = 1;
        events      = handles.events;
        save(fullfile(handles.filedir,handles.filename), 'fileID', 'version', 'events');
    end
end

function SaveCheck(hObject, eventdata, handles)
    if (~handles.isSaved)
        buttonChosen = questdlg('Save changes?', 'MakeMovieGUI', 'Yes', 'No', 'Yes');
        if (strcmp(buttonChosen,'Yes'))
            Save_Callback(hObject, eventdata, handles);
        end
    end
end

function New_Callback(hObject, eventdata, handles)
    SaveCheck(hObject, eventdata, handles);
    handles.events = [];
    set(handles.EventList, 'String', []);
    set(handles.EventList, 'Value', 1);
    handles.filename = [];
    handles.filedir = [];
    handles = UpdateFileStatus(handles, true);
    guidata(hObject, handles);
end

function Open_Callback(hObject, eventdata, handles)
    [filename filedir] = uigetfile('*.mat');
    if (filename == 0), return; end

    file = load(fullfile(filedir,filename));
    if (strcmp(file.fileID, 'EventsList'))
        if (file.version <= 1)
            handles.filename    = filename;
            handles.filedir     = filedir;
            handles.events      = {};
            set(handles.EventList, 'String', []);
            set(handles.EventList, 'Value', 1);

            for i = 1:length(file.events)
                handles = InsertEvent(hObject, handles, file.events{i}, i);
            end

            handles = UpdateFileStatus(handles, true);
            guidata(hObject, handles);
        else
            errordlg('Incompatible eventsList file.');
        end
    else
        errordlg('Corrupt/invalid eventsList file.');
    end
end

function Save_Callback(hObject, eventdata, handles)
    if (isempty(handles.filename))
        SaveAs_Callback(hObject, eventdata, handles);
    else
        handles = UpdateFileStatus(handles, true);
        SaveData(handles);
        guidata(hObject, handles);
    end
end

function SaveAs_Callback(hObject, eventdata, handles)
    [filename, filedir]   = uiputfile('brainMovie.mat');
    if (filename == 0), return; end

    handles.filename    = filename;
    handles.filedir     = filedir;
    handles             = UpdateFileStatus(handles, true);

    SaveData(handles);
    guidata(hObject, handles);
end

function Quit_Callback(hObject, eventdata, handles)
    SaveCheck(hObject, eventdata, handles);
    close(handles.MakeMovieGUI);
end

function EventList_KeyPressFcn(hObject, eventdata, handles)
    if (~isempty(cellfind(eventdata.Modifier, 'control')) || ...
            ~isempty(cellfind(eventdata.Modifier, 'command'))) % done to support mac as well as linux
        switch (eventdata.Key)
            case 'x'
                CutEvent_Callback(hObject, eventdata, handles);
            case 'c'
                CopyEvent_Callback(hObject, eventdata, handles);
            case 'v'
                PasteEvent_Callback(hObject, eventdata, handles);
            case 'm'
                EditEvent_Callback(hObject, eventdata, handles);
            case 'd'
                DeleteEvent_Callback(hObject, eventdata, handles);
            case 'n'
                New_Callback(hObject, eventdata, handles);
            case 'o'
                Open_Callback(hObject, eventdata, handles);
            case 's'
                Save_Callback(hObject, eventdata, handles);
            case 'q'
                Quit_Callback(hObject, eventdata, handles);
            case 'p'
                RenderMovie_Callback(hObject, eventdata, handles, 1);
            case 'e'
                RenderMovie_Callback(hObject, eventdata, handles, 0);
            otherwise
        end
    else
        switch (eventdata.Key)
            case 'a'
                MoveUpButton_Callback(hObject, eventdata, handles);
            case 'z'
                MoveDownButton_Callback(hObject, eventdata, handles);
            case 'delete'
                DeleteEvent_Callback(hObject, eventdata, handles);
            otherwise
        end
    end
end

function CutEvent_Callback(hObject, eventdata, handles)
    if (length(handles.events) < 1), return; end
    selected = get(handles.EventList, 'Value');
    handles.clipboard = handles.events{selected};
    DeleteEvent_Callback(hObject, eventdata, handles);
end

function CopyEvent_Callback(hObject, eventdata, handles)
    if (length(handles.events) < 1), return; end
    selected = get(handles.EventList, 'Value');
    handles.clipboard = handles.events{selected};
    handles = TurnOffInvalidButtons(handles);
    guidata(hObject, handles);
end

function PasteEvent_Callback(hObject, eventdata, handles)
    if (~isempty(handles.clipboard))
        if (length(handles.events) < 1)
            handles = InsertEvent(hObject, handles, handles.clipboard, 1);
        else
            selected = get(handles.EventList, 'Value');
            handles = InsertEvent(hObject, handles, handles.clipboard, selected + 1);
        end
        guidata(hObject, handles);
    end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function rank = cel2rank(cel)
%CEL2RANK Convert temperature from celsius to rankine
%
%  rank = CEL2RANK(cel) converts temperature from celsius to rankine.
%
%  See also RANK2CEL, CEL2FAR, CEL2KEL.

rank = (cel + 273.15) * 1.8;

author: Sable

repo: mcbench-benchmarks


text:
// BLAS 3 - Dense Matrix-Dense Matrix multiplication
// Constants are in GEMM convention; (MxK matrix) x (KxN matrix) = (MxN matrix)
Constant SzM 100;
Constant SzN 100;
Constant SzK 100;

// Optimize for data reuse
// Constant MTileSz 10;
// Constant NTileSz 50;
// Constant KTileSz 100;

// Optimize for throughput
Constant MTileSz 1;
Constant NTileSz 1;
Constant KTileSz 100;

Network BLAS3 {
	Layer BLAS {
		Type: GEMM
		Dimensions { K: SzK, M: SzM, N: SzN }
		Dataflow {
			TemporalMap(NTileSz,NTileSz) N;
			SpatialMap(MTileSz,MTileSz) M;
			TemporalMap(KTileSz,KTileSz) K;
			Cluster(KTileSz, P);
			TemporalMap(NTileSz,NTileSz) N;
			TemporalMap(MTileSz,MTileSz) M;
			SpatialMap(1,1) K;
		}
	}
}

// Network BLAS3 {
//	Layer BLAS {
//		Type: CONV
//		Stride { X: 1, Y: 1 }
//		Dimensions { K: SzN, C: 1, R: 1, S: SzK, Y: SzM, X: SzK }
//		Dataflow {
//			TemporalMap(NTileSz,NTileSz) K;
//			SpatialMap(MTileSz,MTileSz) Y;
//			TemporalMap(SzK,SzK) X;
//			TemporalMap(SzK,SzK) S;
//			Cluster(1,P);
//			SpatialMap(KTileSz,KTileSz) X;
//			TemporalMap(KTileSz,KTileSz) S;
//			Cluster(KTileSz,P);
//			SpatialMap(1,1) X;
//			SpatialMap(1,1) S;
//		}
// 	}
// }

author: maestro-project

repo: maestro


text:
function [VV,FF,SS] = sqrt3(V,F,varargin)
  % SQRT3 Perform sqrt-3 subdivision. After n iterations of loop subivision, the
  % resulting mesh will have
  %   3^n |F| faces
  %   (3^(n+1)-3)/2 |F| + |E| edges
  %   (3^(n+1)-3)/6 |F| + |V| vertices
  %
  % [VV,FF,SS] = sqrt3(V,F)
  %
  % Inputs:
  %   V  #V by dim list of vertex positions
  %   F  #F by simplex-size list of simplex indices
  %   Optional:
  %     'Iterations' followed by number of recursive calls {1}
  % Outputs:
  %   VV  #VV by 3 list of output vertex positions
  %   FF  #FF by 3 list of triangle indices into VV
  %   SS  #VV by #V sparse matrix so that VV = SS*V

  % default values
  iters = 1;
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Iterations'}, ...
    {'iters'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  SS = speye(size(V,1),size(V,1));
  VV = V;
  FF = F;
  for iter = 1:iters
    n = size(SS,1);
    m = size(FF,1);
    C = n+(1:m)';
    E = edges(FF);
    A = adjacency_matrix(FF);
    [AI,AJ] = find(A);
    val = sum(A,1);
    An = (4-2*cos(2*pi./val))./9;
    Andval = An./val;
    A = sparse(AI,AJ,Andval(AI),n,n) + (diag(sparse(1-An)));
    S = [A;sparse(repmat(1:m,3,1)',FF,1/3,m,n)];
    SS = S*SS;
    FF = [FF(:,2) FF(:,3) C;FF(:,3) FF(:,1) C;FF(:,1) FF(:,2) C];
    %BC = barycenter(V,F);
    %VV = [V;BC];
    FF = flip_edges(FF,E,'NonConflicting',true);
  end
  VV = SS*V;

  % F₁ ← 3F₀
  % E₁ ← E₀ + 3F₀
  % V₁ ← V₀ + F₀

end


text:
function [sos,g] = driving_function_imp_nfchoa_pw(N,R,conf)
%DRIVING_FUNCTION_IMP_NFCHOA_PW second-order section representation for a
%plane wave in NFC-HOA
%
%   Usage: sos = driving_function_imp_nfchoa_pw(N,R,conf)
%
%   Input parameters:
%       N       - order of spherical Hankel function
%       R       - radius of secondary source array / m
%       conf    - configuration struct (see SFS_config)
%
%   Output parameters:
%       sos     - second-order section representation
%       g       - scalar gain factor
%
%   See also: sound_field_imp, sound_field_imp_nfchoa,
%       driving_function_imp_nfchoa
%
%   References:
%       Spors, Kuscher, Ahrens (2011) - "Efficient realization of model-based
%       rendering for 2.5-dimensional near-field compensated higher order
%       Ambisonics", IEEE Workshop on Applications of Signal Processing to Audio
%       and Acoustics (WASPAA), pp. 61-64,
%       https://doi.org/10.1109/ASPAA.2011.6082325

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 3;
nargmax = 3;
narginchk(nargmin,nargmax);
isargpositivescalar(N,R);
isargstruct(conf);

%% ===== Configuration ==================================================
c = conf.c;
dimension = conf.dimension;
driving_functions = conf.driving_functions;

%% ===== Computation =====================================================
% Find spherical Hankel function zeros
[z,p] = sphbesselh_zeros(N);

% Get the delay and weighting factors
if strcmp('2D',dimension)

    % === 2-Dimensional ==================================================

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        to_be_implemented;
    otherwise
        error(['%s: %s, this type of driving function is not implemented', ...
            'for a 2D plane wave.'],upper(mfilename),driving_functions);
    end

elseif strcmp('2.5D',dimension) || strcmp('3D',dimension)

    % === 2.5- & 3-Dimensional ==========================================

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        % 2.5D for a plane wave as source model
        %
        [sos, g] = zp2sos(p,z*c/R,1,'down','none');
        g = g * (-1)^abs(N) * 4*pi * R;
        %
        % Compare Spors et al. (2011), eq. (10)
        %
    otherwise
        error(['%s: %s, this type of driving function is not implemented', ...
            'for a 2.5D plane wave.'],upper(mfilename),driving_functions);
    end

else
    error('%s: the dimension %s is unknown.',upper(mfilename),dimension);
end

author: sfstoolbox

repo: sfs-matlab


text:
function [ output_meta ] = meta2sicd_xhdr( xhdr )
%META2SICD_XHDR Convert recognized extended header data to sicd metadata
%structure
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

output_meta = struct();
xhdr_fieldnames = fieldnames(xhdr);
for i=1:length(xhdr_fieldnames)
    if exist(strcat('meta2sicd_',xhdr_fieldnames{i},'.m'),'file');
        eval(['output_meta = setstructfields(output_meta, meta2sicd_'...
            xhdr_fieldnames{i} '(xhdr.' xhdr_fieldnames{i} '));']);
    end
end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function y = vl_nnnormalizelp(x,dzdy,varargin)
%VL_NNNORMALIZELP  CNN Lp normalization
%   Y = VL_NNNORMALIZELP(X) normalizes in Lp norm each column of
%   features in the array X, along its third dimension:
%
%       Y(i,j,k) = X(i,j,k) / sum_q (X(i,j,q).^p + epsilon)^(1/p)
%
%   [Y,N] = VL_NNNORMALIZELP(X) returns the array N containing the
%   computed norms.
%
%   DZDX = VL_NNNORMALIZELP(X, DZDY) computes the derivative of the
%   function with respect to X projected onto DZDY.
%
%   VL_NNNORMALIZE(___, 'opts', val, ...) takes the following options:
%
%   `exponent`:: 2
%      The exponent of the Lp norm. Warning: currently only even
%      exponents are supported.
%
%   `p`:: same as exponent
%
%   `epsilon`:: 0.01
%      The constant added to the sum of p-powers before taking the
%      1/p square root (see the formula above).
%
%   `dimensions`:: [3]
%      The list of dimensions along wich to operate. By default,
%      normalization is along the third dimension, usually
%      corresponding to feature channels.
%
%   `spatial`:: `false`
%      If `true`, sum along the two spatial dimensions instead of
%      along the feature channels. This is the same as setting
%      `dimensions` to [1,2].
%
%   See also: VL_NNNORMALIZE().

opts.epsilon = 1e-2 ;
opts.p = 2 ;
opts.spatial = false ;
opts = vl_argparse(opts, varargin, 'nonrecursive') ;

if ~opts.spatial
  massp = sum(x.^opts.p,3) + opts.epsilon ;
else
  massp = sum(sum(x.^opts.p,1),2) + opts.epsilon ;
end
mass = massp.^(1/opts.p) ;
y = bsxfun(@rdivide, x, mass) ;

if nargin < 2 || isempty(dzdy)
  return ;
else
  dzdy = bsxfun(@rdivide, dzdy, mass) ;
  if ~opts.spatial
    tmp = sum(dzdy .* x, 3) ;
  else
    tmp = sum(sum(dzdy .* x, 1),2);
  end
  y = dzdy - bsxfun(@times, tmp, bsxfun(@rdivide, x.^(opts.p-1), massp)) ;
end

repo: UDT


text:
function  matRad_plotAxisLabels(axesHandle,ct,plane,slice,defaultFontSize,tickdist)
% matRad function to plot x and y labels denoting the ct dimensions
% according to the selected plane
%
% call
%   matRad_plotAxisLabels(axesHandle,ct,plane,slice,defaultFontSize)
%   matRad_plotAxisLabels(axesHandle,ct,plane,slice,defaultFontSize,
%   tickdist)
%
% input
%   axesHandle         handle to axes the slice should be displayed in
%   ct                 matRad ct structure
%   plane              plane view (coronal=1,sagittal=2,axial=3)
%   slice              slice in the selected plane of the 3D cube
%   defaultFontSize    default font size as double value
%
% output
%   -
%
% References
%   -
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright 2017 the matRad development team.
%
% This file is part of the matRad project. It is subject to the license
% terms in the LICENSE file found in the top-level directory of this
% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
% of the matRad project, including this file, may be copied, modified,
% propagated, or distributed except according to the terms contained in the
% LICENSE file.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('defaultFontSize','var') || isempty(defaultFontSize)
    defaultFontSize = 12;
end

if ~exist('tickdist','var') || isempty(tickdist)
    tickdist = 50;
end
%% Set axis labels and plot iso center
if  plane == 3% Axial plane
    if ~isempty(ct.resolution.x) && ~isempty(ct.resolution.y)
        set(axesHandle,'XTick',0:tickdist/ct.resolution.x:1000);
        set(axesHandle,'YTick',0:tickdist/ct.resolution.y:1000);
        set(axesHandle,'XTickLabel',0:tickdist:1000*ct.resolution.x);
        set(axesHandle,'YTickLabel',0:tickdist:1000*ct.resolution.y);
        xlabel(axesHandle,'x [mm]','FontSize',defaultFontSize)
        ylabel(axesHandle,'y [mm]','FontSize',defaultFontSize)
        title(axesHandle,['axial plane z = ' num2str(ct.resolution.z*slice) ' [mm]'],'FontSize',defaultFontSize)
    else
        xlabel(axesHandle,'x [voxels]','FontSize',defaultFontSize)
        ylabel(axesHandle,'y [voxels]','FontSize',defaultFontSize)
        title(axesHandle,'axial plane','FontSize',defaultFontSize)
    end
elseif plane == 2 % Sagittal plane
    if ~isempty(ct.resolution.y) && ~isempty(ct.resolution.z)
        set(axesHandle,'XTick',0:tickdist/ct.resolution.z:1000)
        set(axesHandle,'YTick',0:tickdist/ct.resolution.y:1000)
        set(axesHandle,'XTickLabel',0:tickdist:1000*ct.resolution.z)
        set(axesHandle,'YTickLabel',0:tickdist:1000*ct.resolution.y)
        xlabel(axesHandle,'z [mm]','FontSize',defaultFontSize);
        ylabel(axesHandle,'y [mm]','FontSize',defaultFontSize);
        title(axesHandle,['sagittal plane x = ' num2str(ct.resolution.y*slice) ' [mm]'],'FontSize',defaultFontSize)
    else
        xlabel(axesHandle,'z [voxels]','FontSize',defaultFontSize)
        ylabel(axesHandle,'y [voxels]','FontSize',defaultFontSize)
        title(axesHandle,'sagittal plane','FontSize',defaultFontSize);
    end
elseif plane == 1 % Coronal plane
    if ~isempty(ct.resolution.x) && ~isempty(ct.resolution.z)
        set(axesHandle,'XTick',0:tickdist/ct.resolution.z:1000)
        set(axesHandle,'YTick',0:tickdist/ct.resolution.x:1000)
        set(axesHandle,'XTickLabel',0:tickdist:1000*ct.resolution.z)
        set(axesHandle,'YTickLabel',0:tickdist:1000*ct.resolution.x)
        xlabel(axesHandle,'z [mm]','FontSize',defaultFontSize)
        ylabel(axesHandle,'x [mm]','FontSize',defaultFontSize)
        title(axesHandle,['coronal plane y = ' num2str(ct.resolution.x*slice) ' [mm]'],'FontSize',defaultFontSize)
    else
        xlabel(axesHandle,'z [voxels]','FontSize',defaultFontSize)
        ylabel(axesHandle,'x [voxels]','FontSize',defaultFontSize)
        title(axesHandle,'coronal plane','FontSize',defaultFontSize)
    end
end

end

repo: matRad


text:
% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%%begin

% A serial link manipulator comprises a series of links.  Each link is described
% by four Denavit-Hartenberg parameters.
%
% Let's define a simple 2 link manipulator.  The first link is

L1 = Link('d', 0, 'a', 1, 'alpha', pi/2)

% The Link object we created has a number of properties
L1.a
L1.d

% and we determine that it is a revolute joint
L1.isrevolute

% For a given joint angle, say q=0.2 rad, we can determine the link transform
% matrix
L1.A(0.2)

% The second link is
L2 = Link('d', 0, 'a', 1, 'alpha', 0)

% Now we need to join these into a serial-link robot manipulator

bot = SerialLink([L1 L2], 'name', 'my robot')
% The displayed robot object shows a lot of details.  It also has a number of
% properties such as the number of joints
bot.n

% Given the joint angles q1 = 0.1 and q2 = 0.2 we can determine the pose of the
% robot's end-effector

bot.fkine([0.1 0.2])
% which is referred to as the forward kinematics of the robot.  This, and the
% inverse kinematics are covered in separate demos.

% Finally we can draw a stick figure of our robot

bot.plot([0.1 0.2])

author: petercorke

repo: robotics-toolbox-matlab


text:
function [sigma,mu,var_dirs,var_amts] = scm_ex(Y,A,R,sigma0,mu0,sigma_max)
%SCM_EX Summary of this function goes here
%   Detailed explanation goes here
[N,M] = size(A);
[~,B] = size(R);
sigma_min = 1e-6;

C = kron((A'*A),eye(B));
S = (mu0^2/(sigma0^2))*eye(M*B,M*B);

h = (Y-A*R)'*A;
h = h(:);
gamma = 1/(mu0^2);
lsq = sum(sum((Y-A*R).^2));
gamma0Inv = 1/(N*B)*lsq;
obj_fun_sigma = @(S,gamma) gamma*lsq - gamma*h'*((S+C)\h) ...
    + logdet(S+C) - logdet(S) - N*B*log(gamma);
calc_ignores = @(S,gamma) [(gamma*h'*((S+C)\h)) / (gamma*lsq), ...
    (logdet(S+C) - logdet(S)) / (gamma*lsq)];
disp(['Inital ignores are: ',num2str(calc_ignores(S,gamma))]);
if 0
    % U1'*N1 has a much smaller frobenius norm than that of N1
    global A_gt; global R_gt;
    N = Y-A*R; N1 = Y-A_gt*R_gt;
    [U,D,V] = svd(A,0); [U1,D1,V1] = svd(A_gt,0);
    sum(sum(N.^2)), sum(sum((U'*N).^2))
    sum(sum(N1.^2)), sum(sum((U1'*N1).^2))
end
Es = obj_fun_sigma(S,gamma);
delta_t0 = 1e-9;

T = S + C;
Th = T\h;

disp('Start estimating the uncertainty...');
for i = 1:200
    % update S
    S2_old = S;

    GammaTinv = gamma*Th*Th';
    Tinv = inv(T);
    E0 = obj_fun_sigma(S,gamma);
    delta_t = delta_t0;

    der = zeros(M*B,M*B);
    for j = 1:M
        inds = (j-1)*B+1:j*B;
        Sj = S(inds,inds);
        der(inds,inds) = GammaTinv(inds,inds) - inv(Sj) + Tinv(inds,inds);
    end

    while 1
        S2_new = S - delta_t*der;
        S2_new = (S2_new + S2_new')/2;
        for j = 1:M
            inds = (j-1)*B+1:j*B;
            S3 = S2_new(inds,inds);

            [V,D] = eig(S3);
            d = diag(D);

            epsilon = 1/(gamma*sigma_max^2);
            infinity = 1/(gamma*sigma_min^2);
            d(d<epsilon) = epsilon;
%                 d(d>infinity) = infinity;

            D = diag(d);
            S2_new(inds,inds) = V*D*V';
        end

        E0 = [E0; obj_fun_sigma(S2_new,gamma)];
        if E0(end) >= E0(end-1)
            S = S2_old;
            break;
        else
            S2_old = S2_new;
            delta_t = delta_t*10;
        end
    end
    % Update gamma
    T = S + C;
    Th = T\h;
    gamma = gamma0Inv - 1/(N*B)*h'*Th;
    gamma = 1/gamma;
    % Check the amounts of ignored terms for estimating A and R
    if 1
        total_amt = gamma*lsq;
        amt1 = gamma*h'*((S+C)\h);
        amt2 = logdet(S+C) - logdet(S);
        ratio1(i) = amt1/total_amt;
        ratio2(i) = amt2/total_amt;
    end
    % Check convergence
    Es = [Es;obj_fun_sigma(S,gamma)];
    if Es(end-1) - Es(end) <= 1e-6*(Es(1) - Es(2))
        disp(['Finishes at iteration ',num2str(i),', current energy: ',num2str(Es(end))]);
        disp(['Final ignores are: ',num2str(calc_ignores(S,gamma))]);
        break;
    end

    if mod(i,5) == 0
        disp(['Process iteration ',num2str(i),', current energy: ',num2str(Es(end))]);
    end
end

mu = gamma^(-1/2);
d = (1/mu) * ones(B,1);

[sigma,var_dirs,var_amts] = calc_uncertainty_range(S,d);

author: zhouyuanzxcv

repo: Hyperspectral


text:
function keypointsAll= inf2keypointsImgBP(expidx,parts,annolist_gt)

fprintf('inf2keypointsImgBP()\n');

p = rcnn_exp_params(expidx);
conf = rcnn_config('sub_dir', '/cachedir/test', 'exp_dir', [p.expDir '/' p.shortName]);
if (isfield(p,'infDir'))
    infDir = p.infDir;
else
    infDir = [conf.cache_dir '/inference'];
end

keypointsAll = repmat(struct('imgname','','det',nan(16,3)), length(annolist_gt), 1);

for imgidx = 1:length(annolist_gt)

    fprintf('.');

    fnameInf = [infDir '/imgidx_' padZeros(num2str(imgidx-1),5) '.mat'];
%     if (~exist(fnameInf,'file'))
%         continue;
%     end
    load(fnameInf,'predAll');

    keypointsAll(imgidx).imgname = annolist_gt(imgidx).image.name;
    assert(length(p.pidxs) == length(predAll));

    for i = 1:length(p.pidxs)
        pidx = p.pidxs(i);
        % part is a joint
        assert(parts(pidx+1).pos(1) == parts(pidx+1).pos(2));
        jidx = parts(pidx+1).pos(1);

        det = predAll{i};
%         [val,id] = max(det(:,4)); % for not swapped minus factor graph (expidx 326)
        [val,id] = max(det(:,5));

        x = det(id,1);
        y = det(id,2);
        keypointsAll(imgidx).det(jidx+1,:) = [[x y] val];
    end
    if (~mod(imgidx, 100))
        fprintf(' %d/%d\n',imgidx,length(keypointsAll));
    end
end
fprintf(' done\n');

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function J = warp_jacobian(nx, ny, warp, transform)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%J = WARP_JACOBIAN(NX, NY, WARP, TRANSFORM)
% This function computes the jacobian J of warp transform with respect
% to parameters. In case of homography/euclidean transform, the jacobian depends on
% the parameter values, while in affine/translation case is totally invariant.
%
% Input variables:
% NX:           the x-coordinate values of the horizontal side of ROI (i.e. [xmin:xmax]),
% NY:           the y-coordinate values of vertical side of ROI (i.e. [ymin:ymax]),
% WARP:         the warp transform (used only in homography and euclidean case),
% TRANSFORM:    the type of adopted transform
% {'affine''homography','translation','euclidean'}
%
% Output:
% J:            The jacobian matrix J
%--------------------------------------
% $ Ver: 1.3, 13/5/2012,  released by Georgios D. Evangelidis, INRIA, FRANCE
% Email: georgios.evangelidis@inria.fr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

snx=length(nx);
sny=length(ny);

Jx=nx(ones(1,sny),:);
Jy=ny(ones(1,snx),:)';
J0=0*Jx;
J1=J0+1;

switch lower(transform)
    case 'homography'

    xy=[Jx(:)';Jy(:)';ones(1,snx*sny)];

    %3x3 matrix transformation
    A = warp;
    A(3,3) = 1;

    % new coordinates
    xy_prime = A * xy;

    % division due to homogeneous coordinates
    xy_prime(1,:) = xy_prime(1,:)./xy_prime(3,:);
    xy_prime(2,:) = xy_prime(2,:)./xy_prime(3,:);

    den = xy_prime(3,:)';

    Jx(:) = Jx(:) ./ den;
    Jy(:) = Jy(:) ./ den;
    J1(:) = J1(:) ./ den;

    Jxx_prime = Jx;
    Jxx_prime(:) = Jxx_prime(:) .* xy_prime(1,:)';
    Jyx_prime = Jy;
    Jyx_prime(:) = Jyx_prime(:) .* xy_prime(1,:)';

    Jxy_prime = Jx;
    Jxy_prime(:) = Jxy_prime(:) .* xy_prime(2,:)';
    Jyy_prime = Jy;
    Jyy_prime(:) = Jyy_prime(:) .* xy_prime(2,:)';

    J = [Jx, J0, -Jxx_prime, Jy, J0, - Jyx_prime, J1, J0;...
        J0, Jx, -Jxy_prime, J0, Jy, -Jyy_prime, J0, J1];

        case 'affine'

    J = [Jx, J0, Jy, J0, J1, J0;...
        J0, Jx, J0, Jy, J0, J1];
    case 'translation'

    J = [J1, J0;...
        J0, J1];

    case 'euclidean'

        mycos = warp(1,1);
        mysin = warp(2,1);

        Jx_prime = -mysin*Jx - mycos*Jy;
        Jy_prime =  mycos*Jx - mysin*Jy;

        J = [Jx_prime, J1, J0;...
            Jy_prime, J0, J1];

    otherwise
        error('function WARP_JACOBIAN: Unknown transform!');
end


text:
function kern = kernExpandParamTransformSettings(kern, paramtransformsettings)

% KERNEXPANDPARAMTRANSFORMSETTINGS Expand parameters' transform settings to form a kernel structure.
%
% FORMAT
%
% DESC returns a kernel structure filled with the parameter
% transform settings in the given cell array. This is used as a
% helper function to enable parameters to be optimised in, for
% example, the NETLAB optimisation functions.
%
% ARG kern : the kernel structure in which the parameter
% transformation settings are to be placed.
%
% ARG paramtransformsettings : cell array (nParams x 1) of
% parameter transform settings which are to be placed in the
% kernel structure. Each transformation setting can be of any
% form, depending on what information the corresponding
% transformation needs. For example, a transformation setting
% might be a desired output range like [0 100], or it could
% be some more complicated structure.
%
% RETURN kern : kernel structure with the given parameter
% transform settings in the relevant locations.
%
% SEEALSO : kernExtractParam, scg, conjgrad
%
% COPYRIGHT : Neil D. Lawrence, 2003, 2004, 2005, 2006
%
% COPYRIGHT : Jaakko Peltonen, 2011
%
% COPYRIGHT : Antti Honkela, 2011

% KERN

%fprintf(1,'kernExpandParamTransformSettings step1\n');

% If no transformation settings are provided, do nothing, otherwise
% simply call the kernel-specific function.
if ~isempty(paramtransformsettings),
  % Get the handle of the kernel-specific function and call it
  fhandle = str2func([kern.type 'KernExpandParamTransformSettings']);
  kern = fhandle(kern, paramtransformsettings);
end;

if isfield(kern, 'priors')
  if iscell(kern.priors)
    for k=1:length(kern.priors),
      if isfield(kern.priors{k}, 'isBounded') && ...
	    kern.priors{k}.isBounded,
	l = kern.priors{k}.index;
	if length(l) > 1,
	  error('Unable to set bounds to a prior shared by multiple params')
	end
	kern.priors{k} = priorSetBounds(kern.priors{k}, paramtransformsettings{l});
      end
    end
  else
    for k=1:length(kern.priors),
      if isfield(kern.priors(k), 'isBounded') && ...
	    kern.priors(k).isBounded,
	l = kern.priors(k).index;
	if length(l) > 1,
	  error('Unable to set bounds to a prior shared by multiple params')
	end
	kern.priors(k) = priorSetBounds(kern.priors(k), paramtransformsettings{l});
      end
    end
  end
end

%fprintf(1,'kernExpandParamTransformSettings done\n');

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function idx = tri_find_nonmanifold_vertex(tri, x, scimat)
% TRI_FIND_NONMANIFOLD_VERTEX  Find indices of non-manifold vertices in a
% triangular mesh
%
% IDX = tri_find_nonmanifold_vertex(TRI, RCS)
%
%   This function is a hack so that we can identify which vertices
%   meshcheckrepair() identifies as non-manifold.
%
%   TRI is a 3-column matrix. Each row represents the indices of the three
%   vertices that form a triangle. TRI as a whole represents the closed
%   surface.
%
%   RCS is a 3-column matrix. Each row represents the (row, column, slice)
%   coordinates of a vertex in the mesh.
%
%   IDX is a vector with the indices of the vertices in X that are
%   non-manifold.
%
% IDX = tri_find_nonmanifold_vertex(TRI, X, SCIMAT)
%
%   Alternatively, vertex coordinates X can be provided as real world
%   (x,y,z) Cartesian coordinates. SCIMAT is a struct with the offset,
%   voxel size and orientation information (see "help scimat" for details).
%
% See also: scimat_world2index, scimat_index2world.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2013-2014 University of Oxford
% Version: 0.2.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program 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.
%
% This program 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

% check arguments
narginchk(2, 3);
nargoutchk(0, 1);

% defaults
if (nargin < 3 || isempty(scimat))
    scimat = scimat_im2scimat(zeros(2), [1 1 1], [0 0 0]);
end

% we convert the vertex coordinates to indices to avoid loss of precision
% when meshcheckrepair() writes to and reads from OFF file
rcs = round(scimat_world2index(x, scimat));

% meshcheckrepair() is supposed to be used to repair a non-manifold mesh.
% What it actually does is that it duplicates non-manifold vertices, which
% I'm not sure it's a good repair. Anyway, here we are interested only in
% detecting the non-manifold vertices. We are going to do it in a bit of a
% convoluted way, because meshcheckrepair() just duplicated the vertices,
% but doesn't tell which ones they are

% run meshcheckrepair(). This will duplicate any non-manifold vertices We
% use evalc so that we can supress the output to the screen.
[~, rcs2] = evalc('meshcheckrepair(rcs, tri, ''deep'');');

% round values to avoid finite precision errors
rcs2 = round(rcs2);

% assign one index value to each unique vertex
[~, idx] = unique(rcs2, 'rows');

% delete all vertices except for the duplicates
rcs2(idx, :) = [];

% sometimes a vertex may be repetead more than once. We want a list of
% unique vertices
rcs2 = unique(rcs2, 'rows');

% indices of non-manifold vertices
idx = find(ismember(rcs, rcs2, 'rows'));

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function F = extractColumns(F, colIdx)
%EXTRACTCOLUMNS   Extract columns (or rows) of an array-valued CHEBFUN.
%   G = EXTRACTCOLUMNS(F, COLIDX) extracts the columns specified by the row
%   vector COLIDX from the CHEBFUN F so that G = F(:,COLIDX). COLIDX need not be
%   increasing in order or unique, but it must contain only integers in the
%   range [1, M] where F has M columns.
%
%   If F is a row CHEBFUN, then EXTRACTCOLUMNS(F, ROWIDX) behaves as described
%   above, except that extracts the rows of F so G = F(ROWIDX,:).
%
% See also MAT2CELL.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Number of columns (or rows if f.isTransposed) of f:
numCols = numColumns(F);

% Trivial cases:
if ( isempty(F) )
    return

elseif ( ~isnumeric(colIdx) && strcmp(colIdx, ':') )
    return

elseif ( (numel(colIdx) == numel(numCols)) && all(colIdx == 1:numCols) )
    return

elseif ( max(colIdx) > numCols )
    error('CHEBFUN:CHEBFUN:extractColumns:dimensions', ...
        'Index exceeds CHEBFUN dimensions.')

end

if ( numel(F) > 1 )
    % Quasimatrix case:
    F = F(colIdx);

else
    % Array-valued CHEBFUN case:

    % Extract the columns of the FUNs:
    for k = 1:numel(F.funs)
        F.funs{k} = extractColumns(F.funs{k}, colIdx);
    end

    % Extract the columns from the impulses:
    F.pointValues = F.pointValues(:,colIdx);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function C = mtimes(A,B)
%MTIMES Multiplies two tenmat objects.
%
%  C = MTIMES(A,B) computes the product of A and B. The result is a
%  TENMAT object and can be transformed into a tensor.
%
%  C = MTIMES(A,B) is called for the syntax 'A * B' when A or B is a
%  TENMAT object.
%
%   See also TENMAT.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

% Handle scalar input
if ~isa(B,'tenmat') && numel(B) == 1
    C = A;
    C.data = C.data * B;
    return;
end
if ~isa(A,'tenmat') && numel(A) == 1
    C = B;
    C.data = C.data * A;
    return;
end

% Handle matrix input
if ~isa(A,'tenmat')
    A = tenmat(A,1);
end

if ~isa(B,'tenmat')
    B = tenmat(B,1);
end

% Error check
if size(A,2) ~= size(B,1)
    error(['Size mismatch: Number of columns in A is not equal to' ...
	   ' the number of rows in B']);
end

tsiz = [A.tsize(A.rindices) B.tsize(B.cindices)];

if ~isempty(tsiz)
    C = tenmat;
    C.tsize = tsiz;
    C.rindices = 1:length(A.rindices);
    C.cindices = (1:length(B.cindices)) + length(A.rindices);
    C.data = A.data * B.data;
else
    C = A.data * B.data;
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% This is the main script runner for training, it collects the
% training samples, followed by model training

clear

% define the root name of database
root = '../data_preparation/prepared_data/';

% which scales we're doing
scales = [0.25, 0.35, 0.5];

% the data generation parameters
sigma = 1;
num_samples = 5e5;
ratio_neg = 5;
norm = 1;
normalisation_size = 19;

patch_types = {'reg', 'grad'};

% Should a frontal view be created
frontalView = [1];

% The other views to be used
profileViewInds = [2,3,4];
upDownViewInds = [];

version = 'multi_pie';

% the naming of
wild_loc = 'mpie_';

for s=scales
    Train_all(root, frontalView, profileViewInds, upDownViewInds, s,...
        sigma, version, norm, 'ratio_neg', ratio_neg,...
        'num_samples', num_samples, 'data_loc', wild_loc,...
        'patch_types', patch_types, 'normalisation_size', normalisation_size);
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function d = fmris_read_analyze(file,Z,T);
%
% Read in 1 4D or >1 3D Analyze format Image Volumes
%
% Usage
%
% d=Read_Analyze(file);
% Reads in all image data and attributes into the structure d.
%
% d=Read_Analyze(file,Z,T);
% Reads in chosen planes and frames and attributes into the structure d.
% Z and T are vectors e.g. Z=1:31;T=1:12; or Z=[24 26];T=[1 3 5];
%
% (c) Roger Gunn & John Aston
%
% HISTORY:
% 2004.02.15 Bob Dougherty (RFD): fixed minor bug- d.file_name gets
% clobbered, so we use 'file' instead.
% Get machine format on which the file was written:
machineformat = getmachineformat(file);
if ~isempty(machineformat)
   % Read Header Information
   d = Read_Analyze_Hdr(file,machineformat);

   % try to get file precision if it is unknown:
   if strcmp(d.precision,'Unknown')
      d.precision=getprecision(deblank(file),machineformat,d.dim,d.global);
   end
   % Read Image Data
   if ~isempty(d) & ~strcmp(d.precision,'Unknown')

      if nargin==1 | strcmp(d.precision,'uint1')

         % Read in Whole Analyze Volume
         fid = fopen(file,'r',machineformat);
         if fid > -1
            d.data=d.scale*reshape(fread(fid,prod(d.dim),d.precision),d.dim(1),d.dim(2),d.dim(3),d.dim(4));
            fclose(fid);
         else
            errordlg('Check Image File: Existence, Permissions ?','Read Error');
         end;

         if nargin==3
            if all(Z>0)&all(Z<=d.dim(3))&all(T>0)&all(T<=d.dim(4))
               d.data=d.data(:,:,Z,T);
               d.Z=Z;
               d.T=T;
            else
               errordlg('Incompatible Matrix Identifiers !','Read Error');
            end
         end

      elseif nargin==3
         % Read in Chosen Planes and Frames
         if (T(1)~=0)|(Z(1)~=0)

            if all(Z>0)&all(Z<=d.dim(3))&all(T>0)&all(T<=d.dim(4))

               fid = fopen(d.file_name,'r',machineformat);
               if fid > -1
                  d.data=zeros(d.dim(1),d.dim(2),length(Z),length(T));
                  for t=1:length(T)
                     for z=1:length(Z)
                        status=fseek(fid,d.hdr.byte*((T(t)-1)*prod(d.dim(1:3))+(Z(z)-1)*prod(d.dim(1:2))),'bof');
                        d.data(:,:,z,t)=d.scale*fread(fid,[d.dim(1) d.dim(2)],d.precision);
                     end
                  end
                  d.Z=Z;
                  d.T=T;
                  fclose(fid);
               else
                  errordlg('Check Image File: Existence, Permissions ?','Read Error');
               end;

            else
               errordlg('Incompatible Matrix Identifiers !','Read Error');
            end;

         end;
      else
         errordlg('Unusual Number of Arguments','Read Error');
      end;
   else
      if strcmp(d.precision,'Unknown');
         errordlg('Unknown Data Type (Precision?)','Read Error');
      end
   end;
else
   errordlg('Unknown Machine Format','Read Error');
end
% if there is no slice thickness, set it to 6mm:
if d.vox(3)==0
   d.vox(3)=6;
end
d.file_name = file;
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function machineformat=getmachineformat(file);
% Get machine format by reading the d.hdr.dim(1) attribute and
% making sure it is 1, 2, 3 or 4.
machineformat=[];
for mf='nlbdgcas'
   fid = fopen([file(1:(length(file)-3)) 'hdr'],'r',mf);
   if fid > -1
      fseek(fid,40,'bof');
      if any(fread(fid,1,'int16')==1:4)
         machineformat=mf;
         fclose(fid);
         break
      else
         fclose(fid);
      end
   else
      errordlg('Check Header File: Existence, Permissions ?','Read Error');
      break
   end
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function precision=getprecision(file,machineformat,dim,range);
% Get precision by reading a value from the middle of the .img file and
% making sure it is within the global attribute
precisions=['int8   '
        'int16  '
        'int32  '
        'int64  '
        'uint8  '
        'uint16 '
        'uint32 '
        'uint64 '
        'single '
        'float32'
        'double '
        'float64' ];
nbytes=[1 2 4 8 1 2 4 8 4 4 8 8];
middle_vol=dim(1)*dim(2)*floor(dim(3)/2)+dim(1)*round(dim(2)/2)+round(dim(1)/2);
h=dir(file);
n=ceil(h.bytes/prod(dim));
fid = fopen(file,'r',machineformat);
if fid > -1
   for i=1:size(precisions,1)
      if nbytes(i)==n
         status=fseek(fid,middle_vol*n,'bof');
         if status==0
            precision=deblank(precisions(i,:));
            x=fread(fid,10,precision);
            if all(range(1)<=x) & all(x<=range(2))
               return
            end
         end
      end
   end
end
errordlg('Check Header File: Existence, Permissions ?','Read Error');
precision='Unknown';
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d=Read_Analyze_Hdr(file,machineformat);
% Read Analyze Header information into the structure d
% Adapted from John Ashburners spm_hread.m
fid  = fopen([file(1:(length(file)-3)) 'hdr'],'r',machineformat);
if fid > -1

   % read (struct) header_key
   %---------------------------------------------------------------------------
   fseek(fid,0,'bof');

   d.hdr.sizeof_hdr 		= fread(fid,1,'int32');
   d.hdr.data_type  		= deblank(setstr(fread(fid,10,'char'))');
   d.hdr.db_name    		= deblank(setstr(fread(fid,18,'char'))');
   d.hdr.extents    		= fread(fid,1,'int32');
   d.hdr.session_error   	= fread(fid,1,'int16');
   d.hdr.regular    		= deblank(setstr(fread(fid,1,'char'))');
   d.hdr.hkey_un0    		= deblank(setstr(fread(fid,1,'char'))');

   % read (struct) image_dimension
   %---------------------------------------------------------------------------
   fseek(fid,40,'bof');

   d.hdr.dim    			= fread(fid,8,'int16');
   d.hdr.vox_units    		= deblank(setstr(fread(fid,4,'char'))');
   d.hdr.cal_units    		= deblank(setstr(fread(fid,8,'char'))');
   d.hdr.unused1			= fread(fid,1,'int16');
   d.hdr.datatype			= fread(fid,1,'int16');
   d.hdr.bitpix				= fread(fid,1,'int16');
   d.hdr.dim_un0			= fread(fid,1,'int16');
   d.hdr.pixdim				= fread(fid,8,'float');
   d.hdr.vox_offset			= fread(fid,1,'float');
   d.hdr.funused1			= fread(fid,1,'float');
   d.hdr.funused2			= fread(fid,1,'float');
   d.hdr.funused3			= fread(fid,1,'float');
   d.hdr.cal_max			= fread(fid,1,'float');
   d.hdr.cal_min			= fread(fid,1,'float');
   d.hdr.compressed			= fread(fid,1,'int32');
   d.hdr.verified			= fread(fid,1,'int32');
   d.hdr.glmax				= fread(fid,1,'int32');
   d.hdr.glmin				= fread(fid,1,'int32');

   % read (struct) data_history
   %---------------------------------------------------------------------------
   fseek(fid,148,'bof');

   d.hdr.descrip			= deblank(setstr(fread(fid,80,'char'))');
   d.hdr.aux_file			= deblank(setstr(fread(fid,24,'char'))');
   d.hdr.orient				= fread(fid,1,'char');
   d.hdr.origin				= fread(fid,5,'uint16');
   d.hdr.generated			= deblank(setstr(fread(fid,10,'char'))');
   d.hdr.scannum			= deblank(setstr(fread(fid,10,'char'))');
   d.hdr.patient_id			= deblank(setstr(fread(fid,10,'char'))');
   d.hdr.exp_date			= deblank(setstr(fread(fid,10,'char'))');
   d.hdr.exp_time			= deblank(setstr(fread(fid,10,'char'))');
   d.hdr.hist_un0			= deblank(setstr(fread(fid,3,'char'))');
   d.hdr.views				= fread(fid,1,'int32');
   d.hdr.vols_added			= fread(fid,1,'int32');
   d.hdr.start_field		= fread(fid,1,'int32');
   d.hdr.field_skip			= fread(fid,1,'int32');
   d.hdr.omax				= fread(fid,1,'int32');
   d.hdr.omin				= fread(fid,1,'int32');
   d.hdr.smax				= fread(fid,1,'int32');
   d.hdr.smin				= fread(fid,1,'int32');

   fclose(fid);

   % Put important information in main structure
   %---------------------------------------------------------------------------

   d.dim    	  			= d.hdr.dim(2:5)';
   vox 						= d.hdr.pixdim(2:5)';
   if 	vox(4)==0
      vox(4)=[];
   end
   d.vox       				= vox;
   d.vox_units       		= d.hdr.vox_units;
   d.vox_offset	    		= d.hdr.vox_offset;
   scale     				= d.hdr.funused1;
   d.scale     			  	= ~scale + scale;
   d.global					= [d.hdr.glmin d.hdr.glmax];
   d.calib					= [d.hdr.cal_min d.hdr.cal_max];
   d.calib_units			= d.hdr.cal_units;
   d.origin    				= d.hdr.origin(1:3)';
   d.descrip   				= d.hdr.descrip(1:max(find(d.hdr.descrip)));

   switch d.hdr.datatype
   case 1
      d.precision 	= 'uint1';
      d.hdr.byte 	= 0;
   case 2
      d.precision 	= 'uint8';
      d.hdr.byte 	= 1;
   case 4
      d.precision 	= 'int16';
      d.hdr.byte 	= 2;
   case 8
      d.precision 	= 'int32';
      d.hdr.byte 	= 4;
   case 16
      d.precision 	= 'float';
      d.hdr.byte 	= 4;
   case 64
      d.precision 	= 'double';
      d.hdr.byte 	= 8;
   otherwise
      d.precision 	= 'Unknown';
      d.hdr.byte 	= 0;
   end

else
   d=[];
   errordlg('Check Header File: Existence, Permissions ?','Read Error');
end
return


function [sp, reg2sp] = read_sprep(spfile, reg2spfile)
reg2sp=imread(reg2spfile);
B=textread(spfile, '%d');

sz=B(1:2);
sp=reshape(B(3:end),sz(:)');