// Copyright (C) 2025 Avinam Kalma // 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. // // You should have received a copy of the GNU General Public License // along with this program. If not, see http://www.gnu.org/licenses/. #include DEFUN_DLD (__graycomatrix__, args, , "\ -*- texinfo -*-\n\ @deftypefn {Loadable Function} __graycomatrix__(@var{img}, @var{offset}, @var{num_levels}, @var{symmetric})\n\ Compute the gray-level co-occurrence matrix (GLCM) for an image\n\ \n\ @var{img} is a 2D gray-level image.\n\ \n\ @var{offset} is an M-by-2 array specifying the pixel pair offsets for which the GLCM is calculated.\n\ \n\ @var{num_levels} is a scalar specifying the number of gray levels to use when scaling the input image.\n\ \n\ @var{symmetric} is a logical value indicating whether the GLCM should be symmetric.\n\ \n\ This function is internal and should NOT be called directly. Instead use @code{graycomatrix}.\n\ @end deftypefn\n\ ") { // check number of inputs if (args.length() != 4) print_usage(); Matrix scaled_image = args(0).array_value(); Matrix offset = args(1).array_value(); int num_levels = args(2).int_value(); bool symmetric = args(3).bool_value(); // dimensions octave_idx_type rows = scaled_image.rows(); octave_idx_type cols = scaled_image.columns(); octave_idx_type n_offsets = offset.rows(); // prepare output: 3D array flattened in row-major for octave // Create output GLCM NDArray glcm(dim_vector(num_levels, num_levels, n_offsets)); glcm.fill(0.0); // pointers for speed double* image_ptr = scaled_image.fortran_vec(); // column-major double* glcm_ptr = glcm.fortran_vec(); // iterate over offsets for (octave_idx_type k = 0; k < n_offsets; ++k) { int dy = (int)offset(k, 0); int dx = (int)offset(k, 1); int k_offset = k * num_levels * num_levels; // find valid pixels int c_min = std::max(0, -dx); int c_max = std::min(cols, cols - dx); int r_min = std::max(0, -dy); int r_max = std::min(rows, rows - dy); for (octave_idx_type c = c_min; c < c_max; ++c) { int c2 = c + dx; int crows = c * rows; int crows2 = c2 * rows; int r1 = r_min + crows; int r2 = r_min + dy + crows2; for (octave_idx_type r = r_min; r < r_max; ++r) { int gray_level1 = (int)image_ptr[r1++]; int gray_level2 = (int)image_ptr[r2++]; // linear index in 3D: gray_level1 + gray_level2*num_levels + k*num_levels*num_levels glcm_ptr[gray_level1 + gray_level2 * num_levels + k_offset] += 1.0; } } if (symmetric) { for (octave_idx_type i = 0; i < num_levels; ++i) { for (octave_idx_type j = 0; j <= i; ++j) { glcm(i, j, k) += glcm(j, i, k); glcm(j, i, k) = glcm(i, j, k); } } } } return octave_value(glcm); }