From 067ef6c9245a53181c8904316c249d527c3bcbd7 Mon Sep 17 00:00:00 2001 From: leoca Date: Thu, 25 Sep 2025 15:56:54 -0300 Subject: [PATCH 1/3] adding dither function --- scripts/image/dither.m | 377 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 377 insertions(+) create mode 100644 scripts/image/dither.m diff --git a/scripts/image/dither.m b/scripts/image/dither.m new file mode 100644 index 0000000000..bced084f93 --- /dev/null +++ b/scripts/image/dither.m @@ -0,0 +1,377 @@ +## Copyright (C) 2025 Leonardo Araujo +## +## 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. +## +## You should have received a copy of the GNU General Public License +## along with this program; see the file COPYING. If not, see +## . + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{X} = } dither (@var{RGB}, @var{map}) +## @deftypefnx {Function File} {@var{X} = } dither (@var{RGB}, @var{map}, @var{Qm}, @var{Qe}) +## @deftypefnx {Function File} {@var{BW} = } dither (@var{I}) +## +## @code{@var{X} = dither (@var{RGB},@var{map})} creates an indexed image +## approximation, using the color provided in the colormap, and uses dithering +## to increase apparent color resolution. Floyd-Steinberg error filter is used: +## [ x 7] +## [3 5 1] / 16 +## It used a raster scan and no weight renormalization at boundaries. +## The default values are used: @var{Qm}=5, and @var{Qe}=8. +## +## Inputs: +## @var{RGB} is a m x n x 3 array with values in [0, 1] (double) or [0, 255] (uint8). +## @var{map} is c x 3 matrix holding RGB triplets in [0, 1] (double). +## @var{Qm} is the number of quantization bits per axis for inverse colormap (default: 5). +## @var{Qe} is the number of quantization bits for error diffusion (default: 8, max 16). +## +## Output: +## @var{X} is a m x n indexed image (uint8 if c<=256, else uint16) for the +## colormap @var{map} provided. +## +## Example: +## @example +## X = dither (RGB, map); +## @end example +## +## @code{@var{X} = dither (@var{RGB}, @var{map}, @var{Qm}, @var{Qe})} +## +## @var{Qm} is the number of quantization bits along each color axis for the +## inverse colormap. @var{Qm} determines the resolution of this grid along each +## color axis (R, G, B). @var{Qm} defines the precision of the color space +## discretization used to map input RGB values to those colors available in the +## colormap. @var{Qe} is the number of quantization bits for the color space +## error calculations in the Floyd-Steinberg error diffusion algorithm. +## It controls the precision of the error values that are calculated and +## propagated during dithering. If @var{Qe} < @var{Qm}, the error diffusion +## process may lose precision, therefore dithering cannot be performed, and the +## function returns an undithered indexed image. +## +## @code{@var{BW} = dither (@var{I})} converts the grayscale input image @var{I} +## into binary applying dithering in the process. The output image @var{BW} +## is a black and white image where dithering creates the illusion of shades of +## gray. +## +## Ref [1] Floyd, R. W., and Steinberg, L., An Adaptive Algorithm for Spatial +## Gray Scale, International Symposium Digest of Technical Papers, Society for +## Information Displays, 1975, p. 36. +## Ref [2] Ulichney. R., Digital Halftoning, The MIT Press, 1987. +## +## @seealso{rgb2ind, imapprox} +## @end deftypefn + +function X = dither (RGB, map, Qm = 5, Qe = 8) + if (nargin < 1 || nargin > 4 || nargin == 3) + print_usage; + endif + if ndims (RGB) == 2 + RGB = cat (3, RGB, RGB, RGB); % Duplicate grayscale to RGB + if nargin < 2, + map = [0 0 0; 1 1 1]; % binary (black and white) colormap + Qm = 1; + endif + endif + if ndims (RGB) != 3 || size (RGB, 3) != 3 + error('dither: RGB must be an m x n x 3 array.'); + end + if !ismatrix (map) || size (map, 2) != 3 || min (map(:)) < 0 || max (map(:)) > 1 + error('dither: Colormap must be a c x 3 matrix.'); + endif + if nargin > 2, + if Qm < 1 || Qe < 1 || fix (Qm) != Qm || fix (Qe) != Qe + error ('Qm and Qe must be a positive integers.'); + elseif Qe < Qm + warning ('dither: Qe < Qm, returning undithered image.'); + X = zeros (size (RGB, 1), size (RGB, 2), 'uint16'); + for i = 1:size (RGB, 1) + for j = 1:size (RGB, 2) + X(i, j) = rgb2indLUT (RGB(i, j, :), map, Qm); + endfor + endfor + if size (map, 1) <= 256 + X = uint8 (X); + endif + return; + endif + endif + Qe = min (Qe, 16); % Cap Qe to avoid excessive precision + + % Scale RGB and map to [0, 1] + if isa (RGB, 'uint8') + RGB = double (RGB) / 255; + elseif max (RGB(:)) > 1 + RGB = double (RGB) / 255; + end + RGB = max (0, min (1, RGB)); + if max (map(:)) > 1 + map = double (map) / 255; + end + map = max (0, min (1, map)); + + % Initialize output + [h, w, _] = size (RGB); + X = zeros (h, w, 'uint16'); % Indices (1-based) + + % Floyd-Steinberg weights (raster scan, no renormalization) + FSweights = [0 0 7; 3 5 1] / 16; % Sum = 1 + + % neighbor offsets and weights + offsets = [0 1; 1 -1; 1 0; 1 1]; + weights = [FSweights(1, 3), FSweights(2, 1), FSweights(2, 2), FSweights(2, 3)]; + + % Quantization levels for error (Qe) + n_levels = 2^Qe; + error_scale = n_levels - 1; + + % Process pixels in raster order + for i = 1:h + for j = 1:w + % Get current pixel (with accumulated errors) + pixel = RGB(i, j, :); + pixel = reshape (max (0, min (1, pixel)), 1, 3); % Clamp to [0, 1] + + % Quantize to nearest colormap color + id = rgb2indLUT (pixel, map, Qm); + X(i, j) = id; + + % Compute quantization error + chosen_color = map(id+1, :); + error = pixel - chosen_color; % 1x3 + error = round (error * error_scale) / error_scale; % Quantize to Qe bits + + % Diffuse error to neighboring pixels (no renormalization) + for k = 1:length (weights) + ni = i + offsets(k, 1); + nj = j + offsets(k, 2); + if ni <= h && nj >= 1 && nj <= w + % Extract current pixel value as 1x3 + current = reshape (RGB(ni, nj, :), 1, 3); + % Apply weighted error to each channel + new_value = current + error * weights(k); + RGB(ni, nj, :) = reshape (new_value, 1, 3); + endif + endfor + endfor + endfor + + % Convert output to uint8 if colormap size allows + if size (map, 1) <= 256 + X = uint8 (X); + endif +endfunction + +function id = rgb2indLUT (pixel, map, Qm = 5) + % RGB2INDLUT Map an RGB pixel to the nearest colormap index using a lookup table. + % id = RGB2INDLUT (pixel, map) returns the 1-based index of the closest color + % in the colormap 'map' for the input RGB pixel (1x3 vector), using a quantized + % inverse colormap with 2^5 bins per RGB axis. + % id = RGB2INDLUT (pixel, map, Qm) uses Qm bits for quantization per RGB axis. + % + % Inputs: + % pixel: 1x3 vector [R, G, B], values in [0, 1] (double) or [0, 255] (uint8). + % map: c-by-3 matrix, each row an RGB triplet in [0, 1] (double). + % Qm: Number of quantization bits per axis (default: 5). + % + % Output: + % id: Index (1-based) into the colormap 'map' for the closest color. + % + % Notes: + % - Uses a persistent lookup table (LUT) for speed. + % - LUT is recomputed if map or Qm changes. + % - Warns if Qm is too large (>8) due to memory constraints. + % - Assumes input pixel and map are properly scaled (pixel auto-scaled if needed). + + % Validate inputs + if nargin < 2 + error ('rgb2indLUT: Not enough input arguments. Pixel and colormap required.'); + endif + if length (pixel) != 3 + error ('rgb2indLUT: Pixel must be a 1x3 RGB vector.'); + if !isvector (pixel) + [s, i] = sort (size (pixel),'descend'); + pixel = permute (pixel, i); + endif + endif + if !ismatrix (map) || size (map, 2) != 3 + error ('rgb2indLUT: Colormap must be a c-by-3 matrix.'); + endif + if nargin < 3 + Qm = 5; % Default quantization bits + end + if !isscalar (Qm) || Qm < 1 || floor (Qm) != Qm + error ('rgb2indLUT: Qm must be a positive integer.'); + end + if Qm > 8 + warning ('rgb2indLUT: Qm > 8 may use excessive memory (%d^3 bins).', 2^Qm); + endif + + % Scale pixel to [0, 1] + if isa (pixel, 'uint8') + pixel = double (pixel) / 255; % Convert to [0, 1] + elseif max (pixel(:)) > 1 + pixel = double (pixel) / 255; % Assume [0, 255] if values exceed 1 + endif + pixel = max (0, min (1, pixel)); % Clamp to [0, 1] + + % Ensure map is in [0, 1] + if max (map(:)) > 1 + map = double(map) / 255; + endif + map = max (0, min (1, map)); + + % Persistent variables for LUT + persistent lut last_map last_Qm; + + % Check if we need to recompute the LUT + recompute = isempty (lut) || Qm != last_Qm || !isequal (map, last_map); + + % Number of bins per axis + n_bins = 2^Qm; + + % Scale pixel to [0, n_bins-1] for indexing + bin_idx = round (pixel * (n_bins - 1)) + 1; + + if recompute + % Initialize LUT: n_bins x n_bins x n_bins array of colormap indices + lut = zeros (n_bins, n_bins, n_bins, 'uint16'); + + % Compute bin centers for distance calculations + bin_centers = (0:(n_bins-1))' / (n_bins-1); % [0, 1] range, column vector + [R, G, B] = ndims_grid(n_bins, n_bins, n_bins); % Meshgrid for bin indices + bin_rgb = [bin_centers(R(:)), bin_centers(G(:)), bin_centers(B(:))]; % n_bins^3 x 3 + + % Compute Euclidean distances from each bin to each colormap color + c = size (map, 1); % Number of colors + distances = zeros (n_bins^3, c); + for i = 1:c + diff = bin_rgb - map(i, :); % n_bins^3 x 3 + distances(:, i) = sqrt (sum (diff.^2, 2)); % n_bins^3 x 1 + endfor + + % Find the nearest colormap index (1-based) for each bin + [_, indices] = min (distances, [], 2); + lut(:) = indices; % Assign to LUT + + % Update cached parameters + last_map = map; + last_Qm = Qm; + endif + + % Look up the colormap index + id = lut(bin_idx(1), bin_idx(2), bin_idx(3)) - 1; +endfunction + +function [X, Y, Z] = ndims_grid (nx, ny, nz) + % NDIMS_GRID Create 3D grid indices (emulates meshgrid for 3D). + [x, y, z] = ind2sub ([nx, ny, nz], 1:(nx*ny*nz)); + X = reshape (x, nx, ny, nz); + Y = reshape (y, nx, ny, nz); + Z = reshape (z, nx, ny, nz); +endfunction + +%!demo +%! ## Solid gray +%! +%! I = ones (256)/2; +%! X = dither (I); +%! figure; +%! subplot (121); imshow (I); title ('original'); +%! subplot (122); imshow (double(X)); title ('dithered'); + +%!demo +%! ## Four solid gray levels +%! +%! I = [ones(256,64)/4, ones(256,64)/2, ones(256,64)*3/4, ones(256,64)*7/8]; +%! X = dither (I); +%! figure; +%! subplot (121); imshow (I); title ('original'); +%! subplot (122); imshow (double(X)); title ('dithered'); + +%!demo +%! ## Black-White Gradient +%! +%! I = repmat ([0:255]./255,256,1); +%! X = dither (I); +%! figure; +%! subplot (121); imshow (I); title ('original'); +%! subplot (122); imshow (double(X)); title ('dithered'); + +%!demo +%! ## Color Gradient +%! +%! width = 256; height = 256; +%! upperleft = [1, 0, 0]; % Red +%! upperright = [0, 1, 0]; % Green +%! lowerleft = [0, 0, 1]; % Blue +%! lowerright = [0, 0, 0]; % Black +%! +%! % Create a grid for interpolation +%! [x, y] = meshgrid(linspace(0, 1, width), linspace(0, 1, height)); +%! % Initialize the 3D array for the image +%! image = zeros(height, width, 3); +%! % Calculate the interpolated colors for each point +%! % The logic is a bilinear interpolation of the four corner colors +%! % The first dimension of the `image` matrix is the height (y-axis) and the second is the width (x-axis) +%! image(:, :, 1) = (1 - x) .* (1 - y) * lowerleft(1) + x .* (1 - y) * lowerright(1) + (1 - x) .* y * upperleft(1) + x .* y * upperright(1); +%! image(:, :, 2) = (1 - x) .* (1 - y) * lowerleft(2) + x .* (1 - y) * lowerright(2) + (1 - x) .* y * upperleft(2) + x .* y * upperright(2); +%! image(:, :, 3) = (1 - x) .* (1 - y) * lowerleft(3) + x .* (1 - y) * lowerright(3) + (1 - x) .* y * upperleft(3) + x .* y * upperright(3); +%! +%! % Use the corner colors to define the colormap +%! map = [upperleft; upperright; lowerleft; lowerright]; +%! % Apply dither +%! X = dither (image, map); +%! +%! % Display the results +%! figure; +%! subplot (121); imshow (image); title ('original'); +%! subplot (122); imshow (reshape(map(X(:)+1,:), [size(X) 3])); title ('dithered'); + +%!demo +%! # Lenna +%! url = 'https://upload.wikimedia.org/wikipedia/en/7/7d/Lenna_%28test_image%29.png'; +%! rgb_image = imread(url); +%! map = [226 143 122; 199 127 124; 175 71 82; 230 191 168; 210 100 98; 132 50 81; 94 24 65; 149 97 139] / 255; +%! X = dither (rgb_image, map); +%! I = reshape(map(X(:)+1,:), [size(X) 3]); +%! figure; +%! subplot (121); imshow (rgb_image); title ('original'); +%! subplot (122); imshow (I); title ('dithered'); + +## Test input validation +%!error dither () +%!error dither (permute (1:3,[1 3 2])) +%!error dither (1, 1) +%!error dither (1, 1:3) +%!error dither (1, [0 0 0]') +%!error dither (1, [0 0 0], 0) +%!error dither (1, [0 0 0], 0, 0) +%!error dither (1, [0 0 0], -1, 1) +%!error dither (1, [0 0 0], 1, -1) + +%!test +%! X = dither (0, [0 0 0; 1 1 1], 1, 1); +%! assert (X, uint8(0)) + +%!test +%! X = dither (1, [0 0 0; 1 1 1], 1, 1); +%! assert (X, uint8(1)) + +%!test +%! X = dither (repmat(ones(3)/2,1,1,3), [0 0 0; 1 1 1], 4, 4); +%! assert (X, uint8([1 0 1; 0 1 0; 1 0 1])) + +%!test +%! X = dither (repmat(ones(3)/4,1,1,3), [0 0 0; 1 1 1], 4, 4); +%! assert (X, uint8([0 0 0; 0 1 0; 0 0 0])) + +%!test +%! X = dither (repmat(ones(3)*3/4,1,1,3), [0 0 0; 1 1 1], 4, 4); +%! assert (X, uint8([1 1 1; 1 0 1; 1 1 1])) From 32b1785389360b84e320e6d680694933d91b1c3e Mon Sep 17 00:00:00 2001 From: Leonardo Araujo Date: Tue, 14 Oct 2025 14:52:43 -0300 Subject: [PATCH 2/3] Update dither.m --- scripts/image/dither.m | 34 ++++++++++++++++++---------------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/scripts/image/dither.m b/scripts/image/dither.m index bced084f93..aba50c307d 100644 --- a/scripts/image/dither.m +++ b/scripts/image/dither.m @@ -18,23 +18,25 @@ ## @deftypefn {Function File} {@var{X} = } dither (@var{RGB}, @var{map}) ## @deftypefnx {Function File} {@var{X} = } dither (@var{RGB}, @var{map}, @var{Qm}, @var{Qe}) ## @deftypefnx {Function File} {@var{BW} = } dither (@var{I}) +## Quantize an image, using dithering to increase the apparent color resolution. ## ## @code{@var{X} = dither (@var{RGB},@var{map})} creates an indexed image -## approximation, using the color provided in the colormap, and uses dithering -## to increase apparent color resolution. Floyd-Steinberg error filter is used: +## approximation. It uses the color provided in the colormap, and uses dithering +## to increase apparent color resolution. Floyd-Steinberg error filter is: ## [ x 7] ## [3 5 1] / 16 -## It used a raster scan and no weight renormalization at boundaries. +## It uses a raster scan and no weight renormalization at boundaries. ## The default values are used: @var{Qm}=5, and @var{Qe}=8. ## -## Inputs: -## @var{RGB} is a m x n x 3 array with values in [0, 1] (double) or [0, 255] (uint8). -## @var{map} is c x 3 matrix holding RGB triplets in [0, 1] (double). +## @var{RGB} is a mxnx3 array with values in [0, 1] (double) or [0, 255] (uint8). +## +## @var{map} is cx3 matrix holding RGB triplets in [0, 1] (double). +## ## @var{Qm} is the number of quantization bits per axis for inverse colormap (default: 5). +## ## @var{Qe} is the number of quantization bits for error diffusion (default: 8, max 16). ## -## Output: -## @var{X} is a m x n indexed image (uint8 if c<=256, else uint16) for the +## @var{X} is a mxn indexed image (uint8 if c<=256, else uint16) for the ## colormap @var{map} provided. ## ## Example: @@ -52,7 +54,7 @@ ## error calculations in the Floyd-Steinberg error diffusion algorithm. ## It controls the precision of the error values that are calculated and ## propagated during dithering. If @var{Qe} < @var{Qm}, the error diffusion -## process may lose precision, therefore dithering cannot be performed, and the +## process may lose precision. Therefore dithering cannot be performed, and the ## function returns an undithered indexed image. ## ## @code{@var{BW} = dither (@var{I})} converts the grayscale input image @var{I} @@ -73,20 +75,20 @@ print_usage; endif if ndims (RGB) == 2 - RGB = cat (3, RGB, RGB, RGB); % Duplicate grayscale to RGB + RGB = cat (3, RGB, RGB, RGB); # Duplicate grayscale to RGB if nargin < 2, map = [0 0 0; 1 1 1]; % binary (black and white) colormap Qm = 1; endif endif - if ndims (RGB) != 3 || size (RGB, 3) != 3 - error('dither: RGB must be an m x n x 3 array.'); + if (ndims (RGB) != 3 || size (RGB, 3) != 3) + error ('dither: RGB must be an m x n x 3 array.'); end - if !ismatrix (map) || size (map, 2) != 3 || min (map(:)) < 0 || max (map(:)) > 1 - error('dither: Colormap must be a c x 3 matrix.'); + if (! ismatrix (map) || size (map, 2) != 3 || min (map(:)) < 0 || max (map(:)) > 1) + error ('dither: Colormap must be a c x 3 matrix.'); endif if nargin > 2, - if Qm < 1 || Qe < 1 || fix (Qm) != Qm || fix (Qe) != Qe + if (Qm < 1 || Qe < 1 || fix (Qm) != Qm || fix (Qe) != Qe) error ('Qm and Qe must be a positive integers.'); elseif Qe < Qm warning ('dither: Qe < Qm, returning undithered image.'); @@ -104,7 +106,7 @@ endif Qe = min (Qe, 16); % Cap Qe to avoid excessive precision - % Scale RGB and map to [0, 1] + ## Scale RGB and map to [0, 1] if isa (RGB, 'uint8') RGB = double (RGB) / 255; elseif max (RGB(:)) > 1 From 8741e8af354f63e9e86dccc3d34764ca4781fa3f Mon Sep 17 00:00:00 2001 From: Leonardo Araujo Date: Fri, 17 Oct 2025 13:33:58 -0300 Subject: [PATCH 3/3] Refactor comments and formatting in dither.m Updated comments and formatting in the dither.m file to comply with style guidelines. --- scripts/image/dither.m | 218 ++++++++++++++++++++--------------------- 1 file changed, 109 insertions(+), 109 deletions(-) diff --git a/scripts/image/dither.m b/scripts/image/dither.m index aba50c307d..1052efc3ab 100644 --- a/scripts/image/dither.m +++ b/scripts/image/dither.m @@ -11,7 +11,7 @@ ## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License -## along with this program; see the file COPYING. If not, see +## along with this program; see the file COPYING. If not, see ## . ## -*- texinfo -*- @@ -21,8 +21,8 @@ ## Quantize an image, using dithering to increase the apparent color resolution. ## ## @code{@var{X} = dither (@var{RGB},@var{map})} creates an indexed image -## approximation. It uses the color provided in the colormap, and uses dithering -## to increase apparent color resolution. Floyd-Steinberg error filter is: +## approximation. It uses the color provided in the colormap, and uses dithering +## to increase apparent color resolution. Floyd-Steinberg error filter is: ## [ x 7] ## [3 5 1] / 16 ## It uses a raster scan and no weight renormalization at boundaries. @@ -54,7 +54,7 @@ ## error calculations in the Floyd-Steinberg error diffusion algorithm. ## It controls the precision of the error values that are calculated and ## propagated during dithering. If @var{Qe} < @var{Qm}, the error diffusion -## process may lose precision. Therefore dithering cannot be performed, and the +## process may lose precision. Therefore dithering cannot be performed, and the ## function returns an undithered indexed image. ## ## @code{@var{BW} = dither (@var{I})} converts the grayscale input image @var{I} @@ -74,10 +74,10 @@ if (nargin < 1 || nargin > 4 || nargin == 3) print_usage; endif - if ndims (RGB) == 2 + if (ndims (RGB) == 2) RGB = cat (3, RGB, RGB, RGB); # Duplicate grayscale to RGB if nargin < 2, - map = [0 0 0; 1 1 1]; % binary (black and white) colormap + map = [0 0 0; 1 1 1]; # binary (black and white) colormap Qm = 1; endif endif @@ -87,7 +87,7 @@ if (! ismatrix (map) || size (map, 2) != 3 || min (map(:)) < 0 || max (map(:)) > 1) error ('dither: Colormap must be a c x 3 matrix.'); endif - if nargin > 2, + if (nargin > 2) if (Qm < 1 || Qe < 1 || fix (Qm) != Qm || fix (Qe) != Qe) error ('Qm and Qe must be a positive integers.'); elseif Qe < Qm @@ -98,181 +98,181 @@ X(i, j) = rgb2indLUT (RGB(i, j, :), map, Qm); endfor endfor - if size (map, 1) <= 256 + if (size (map, 1) <= 256) X = uint8 (X); endif return; endif endif - Qe = min (Qe, 16); % Cap Qe to avoid excessive precision + Qe = min (Qe, 16); # Cap Qe to avoid excessive precision ## Scale RGB and map to [0, 1] - if isa (RGB, 'uint8') + if (isa (RGB, 'uint8')) RGB = double (RGB) / 255; - elseif max (RGB(:)) > 1 + elseif (max (RGB(:)) > 1) RGB = double (RGB) / 255; end RGB = max (0, min (1, RGB)); - if max (map(:)) > 1 + if (max (map(:)) > 1) map = double (map) / 255; end map = max (0, min (1, map)); - % Initialize output + ## Initialize output [h, w, _] = size (RGB); - X = zeros (h, w, 'uint16'); % Indices (1-based) + X = zeros (h, w, 'uint16'); # Indices (1-based) - % Floyd-Steinberg weights (raster scan, no renormalization) - FSweights = [0 0 7; 3 5 1] / 16; % Sum = 1 + ## Floyd-Steinberg weights (raster scan, no renormalization) + FSweights = [0 0 7; 3 5 1] / 16; # Sum = 1 - % neighbor offsets and weights + ## neighbor offsets and weights offsets = [0 1; 1 -1; 1 0; 1 1]; weights = [FSweights(1, 3), FSweights(2, 1), FSweights(2, 2), FSweights(2, 3)]; - % Quantization levels for error (Qe) + ## Quantization levels for error (Qe) n_levels = 2^Qe; error_scale = n_levels - 1; - % Process pixels in raster order + ## Process pixels in raster order for i = 1:h for j = 1:w - % Get current pixel (with accumulated errors) + ## Get current pixel (with accumulated errors) pixel = RGB(i, j, :); - pixel = reshape (max (0, min (1, pixel)), 1, 3); % Clamp to [0, 1] + pixel = reshape (max (0, min (1, pixel)), 1, 3); # Clamp to [0, 1] - % Quantize to nearest colormap color + ## Quantize to nearest colormap color id = rgb2indLUT (pixel, map, Qm); X(i, j) = id; - % Compute quantization error + ## Compute quantization error chosen_color = map(id+1, :); - error = pixel - chosen_color; % 1x3 - error = round (error * error_scale) / error_scale; % Quantize to Qe bits + error = pixel - chosen_color; # 1x3 + error = round (error * error_scale) / error_scale; # Quantize to Qe bits - % Diffuse error to neighboring pixels (no renormalization) + ## Diffuse error to neighboring pixels (no renormalization) for k = 1:length (weights) ni = i + offsets(k, 1); nj = j + offsets(k, 2); - if ni <= h && nj >= 1 && nj <= w - % Extract current pixel value as 1x3 + if (ni <= h && nj >= 1 && nj <= w) + ## Extract current pixel value as 1x3 current = reshape (RGB(ni, nj, :), 1, 3); - % Apply weighted error to each channel + ## Apply weighted error to each channel new_value = current + error * weights(k); - RGB(ni, nj, :) = reshape (new_value, 1, 3); + RGB(ni, nj, :) = reshape (new_value, 1, 3); endif endfor endfor endfor - % Convert output to uint8 if colormap size allows + ## Convert output to uint8 if colormap size allows if size (map, 1) <= 256 X = uint8 (X); endif endfunction function id = rgb2indLUT (pixel, map, Qm = 5) - % RGB2INDLUT Map an RGB pixel to the nearest colormap index using a lookup table. - % id = RGB2INDLUT (pixel, map) returns the 1-based index of the closest color - % in the colormap 'map' for the input RGB pixel (1x3 vector), using a quantized - % inverse colormap with 2^5 bins per RGB axis. - % id = RGB2INDLUT (pixel, map, Qm) uses Qm bits for quantization per RGB axis. - % - % Inputs: - % pixel: 1x3 vector [R, G, B], values in [0, 1] (double) or [0, 255] (uint8). - % map: c-by-3 matrix, each row an RGB triplet in [0, 1] (double). - % Qm: Number of quantization bits per axis (default: 5). - % - % Output: - % id: Index (1-based) into the colormap 'map' for the closest color. - % - % Notes: - % - Uses a persistent lookup table (LUT) for speed. - % - LUT is recomputed if map or Qm changes. - % - Warns if Qm is too large (>8) due to memory constraints. - % - Assumes input pixel and map are properly scaled (pixel auto-scaled if needed). - - % Validate inputs - if nargin < 2 - error ('rgb2indLUT: Not enough input arguments. Pixel and colormap required.'); + ## RGB2INDLUT Map an RGB pixel to the nearest colormap index using a lookup table. + ## id = RGB2INDLUT (pixel, map) returns the 1-based index of the closest color + ## in the colormap 'map' for the input RGB pixel (1x3 vector), using a quantized + ## inverse colormap with 2^5 bins per RGB axis. + ## id = RGB2INDLUT (pixel, map, Qm) uses Qm bits for quantization per RGB axis. + ## + ## Inputs: + ## pixel: 1x3 vector [R, G, B], values in [0, 1] (double) or [0, 255] (uint8). + ## map: c-by-3 matrix, each row an RGB triplet in [0, 1] (double). + ## Qm: Number of quantization bits per axis (default: 5). + ## + ## Output: + ## id: Index (1-based) into the colormap 'map' for the closest color. + ## + ## Notes: + ## - Uses a persistent lookup table (LUT) for speed. + ## - LUT is recomputed if map or Qm changes. + ## - Warns if Qm is too large (>8) due to memory constraints. + ## - Assumes input pixel and map are properly scaled (pixel auto-scaled if needed). + + ## Validate inputs + if (nargin < 2) + error ('rgb2indLUT: Not enough input arguments. Pixel and colormap required.'); endif - if length (pixel) != 3 + if (length (pixel) != 3) error ('rgb2indLUT: Pixel must be a 1x3 RGB vector.'); - if !isvector (pixel) + if (! isvector (pixel)) [s, i] = sort (size (pixel),'descend'); pixel = permute (pixel, i); endif endif - if !ismatrix (map) || size (map, 2) != 3 + if (! ismatrix (map) || size (map, 2) != 3) error ('rgb2indLUT: Colormap must be a c-by-3 matrix.'); endif - if nargin < 3 + if (nargin < 3) Qm = 5; % Default quantization bits end - if !isscalar (Qm) || Qm < 1 || floor (Qm) != Qm + if (! isscalar (Qm) || Qm < 1 || floor (Qm) != Qm) error ('rgb2indLUT: Qm must be a positive integer.'); end - if Qm > 8 + if (Qm > 8) warning ('rgb2indLUT: Qm > 8 may use excessive memory (%d^3 bins).', 2^Qm); endif - % Scale pixel to [0, 1] - if isa (pixel, 'uint8') + ## Scale pixel to [0, 1] + if (isa (pixel, 'uint8')) pixel = double (pixel) / 255; % Convert to [0, 1] - elseif max (pixel(:)) > 1 + elseif (max (pixel(:)) > 1) pixel = double (pixel) / 255; % Assume [0, 255] if values exceed 1 endif pixel = max (0, min (1, pixel)); % Clamp to [0, 1] - % Ensure map is in [0, 1] - if max (map(:)) > 1 + ## Ensure map is in [0, 1] + if (max (map(:)) > 1) map = double(map) / 255; endif map = max (0, min (1, map)); - % Persistent variables for LUT + ## Persistent variables for LUT persistent lut last_map last_Qm; - % Check if we need to recompute the LUT - recompute = isempty (lut) || Qm != last_Qm || !isequal (map, last_map); + ## Check if we need to recompute the LUT + recompute = (isempty (lut) || Qm != last_Qm || ! isequal (map, last_map)); - % Number of bins per axis + ## Number of bins per axis n_bins = 2^Qm; - % Scale pixel to [0, n_bins-1] for indexing + ## Scale pixel to [0, n_bins-1] for indexing bin_idx = round (pixel * (n_bins - 1)) + 1; if recompute - % Initialize LUT: n_bins x n_bins x n_bins array of colormap indices + ## Initialize LUT: n_bins x n_bins x n_bins array of colormap indices lut = zeros (n_bins, n_bins, n_bins, 'uint16'); - % Compute bin centers for distance calculations - bin_centers = (0:(n_bins-1))' / (n_bins-1); % [0, 1] range, column vector - [R, G, B] = ndims_grid(n_bins, n_bins, n_bins); % Meshgrid for bin indices - bin_rgb = [bin_centers(R(:)), bin_centers(G(:)), bin_centers(B(:))]; % n_bins^3 x 3 + ## Compute bin centers for distance calculations + bin_centers = (0:(n_bins-1))' / (n_bins-1); # [0, 1] range, column vector + [R, G, B] = ndims_grid (n_bins, n_bins, n_bins); # Meshgrid for bin indices + bin_rgb = [bin_centers(R(:)), bin_centers(G(:)), bin_centers(B(:))]; # n_bins^3 x 3 - % Compute Euclidean distances from each bin to each colormap color - c = size (map, 1); % Number of colors + ## Compute Euclidean distances from each bin to each colormap color + c = size (map, 1); # Number of colors distances = zeros (n_bins^3, c); for i = 1:c - diff = bin_rgb - map(i, :); % n_bins^3 x 3 - distances(:, i) = sqrt (sum (diff.^2, 2)); % n_bins^3 x 1 + diff = bin_rgb - map(i, :); # n_bins^3 x 3 + distances(:, i) = sqrt (sum (diff.^2, 2)); # n_bins^3 x 1 endfor - % Find the nearest colormap index (1-based) for each bin + ## Find the nearest colormap index (1-based) for each bin [_, indices] = min (distances, [], 2); - lut(:) = indices; % Assign to LUT + lut(:) = indices; # Assign to LUT - % Update cached parameters + ## Update cached parameters last_map = map; last_Qm = Qm; endif - % Look up the colormap index + ## Look up the colormap index id = lut(bin_idx(1), bin_idx(2), bin_idx(3)) - 1; endfunction function [X, Y, Z] = ndims_grid (nx, ny, nz) - % NDIMS_GRID Create 3D grid indices (emulates meshgrid for 3D). + ## NDIMS_GRID Create 3D grid indices (emulates meshgrid for 3D). [x, y, z] = ind2sub ([nx, ny, nz], 1:(nx*ny*nz)); X = reshape (x, nx, ny, nz); Y = reshape (y, nx, ny, nz); @@ -295,7 +295,7 @@ %! X = dither (I); %! figure; %! subplot (121); imshow (I); title ('original'); -%! subplot (122); imshow (double(X)); title ('dithered'); +%! subplot (122); imshow (double (X)); title ('dithered'); %!demo %! ## Black-White Gradient @@ -304,37 +304,37 @@ %! X = dither (I); %! figure; %! subplot (121); imshow (I); title ('original'); -%! subplot (122); imshow (double(X)); title ('dithered'); +%! subplot (122); imshow (double (X)); title ('dithered'); %!demo %! ## Color Gradient %! %! width = 256; height = 256; -%! upperleft = [1, 0, 0]; % Red -%! upperright = [0, 1, 0]; % Green -%! lowerleft = [0, 0, 1]; % Blue -%! lowerright = [0, 0, 0]; % Black +%! upperleft = [1, 0, 0]; # Red +%! upperright = [0, 1, 0]; # Green +%! lowerleft = [0, 0, 1]; # Blue +%! lowerright = [0, 0, 0]; # Black %! -%! % Create a grid for interpolation -%! [x, y] = meshgrid(linspace(0, 1, width), linspace(0, 1, height)); -%! % Initialize the 3D array for the image -%! image = zeros(height, width, 3); -%! % Calculate the interpolated colors for each point -%! % The logic is a bilinear interpolation of the four corner colors -%! % The first dimension of the `image` matrix is the height (y-axis) and the second is the width (x-axis) +%! ## Create a grid for interpolation +%! [x, y] = meshgrid (linspace (0, 1, width), linspace (0, 1, height)); +%! ## Initialize the 3D array for the image +%! image = zeros (height, width, 3); +%! ## Calculate the interpolated colors for each point +%! ## The logic is a bilinear interpolation of the four corner colors +%! ## The first dimension of the `image` matrix is the height (y-axis) and the second is the width (x-axis) %! image(:, :, 1) = (1 - x) .* (1 - y) * lowerleft(1) + x .* (1 - y) * lowerright(1) + (1 - x) .* y * upperleft(1) + x .* y * upperright(1); %! image(:, :, 2) = (1 - x) .* (1 - y) * lowerleft(2) + x .* (1 - y) * lowerright(2) + (1 - x) .* y * upperleft(2) + x .* y * upperright(2); %! image(:, :, 3) = (1 - x) .* (1 - y) * lowerleft(3) + x .* (1 - y) * lowerright(3) + (1 - x) .* y * upperleft(3) + x .* y * upperright(3); %! -%! % Use the corner colors to define the colormap +%! ## Use the corner colors to define the colormap %! map = [upperleft; upperright; lowerleft; lowerright]; %! % Apply dither %! X = dither (image, map); %! -%! % Display the results +%! ## Display the results %! figure; %! subplot (121); imshow (image); title ('original'); -%! subplot (122); imshow (reshape(map(X(:)+1,:), [size(X) 3])); title ('dithered'); +%! subplot (122); imshow (reshape (map(X(:)+1,:), [size(X) 3])); title ('dithered'); %!demo %! # Lenna @@ -342,7 +342,7 @@ %! rgb_image = imread(url); %! map = [226 143 122; 199 127 124; 175 71 82; 230 191 168; 210 100 98; 132 50 81; 94 24 65; 149 97 139] / 255; %! X = dither (rgb_image, map); -%! I = reshape(map(X(:)+1,:), [size(X) 3]); +%! I = reshape (map(X(:)+1,:), [size(X) 3]); %! figure; %! subplot (121); imshow (rgb_image); title ('original'); %! subplot (122); imshow (I); title ('dithered'); @@ -360,20 +360,20 @@ %!test %! X = dither (0, [0 0 0; 1 1 1], 1, 1); -%! assert (X, uint8(0)) +%! assert (X, uint8 (0)) %!test %! X = dither (1, [0 0 0; 1 1 1], 1, 1); -%! assert (X, uint8(1)) +%! assert (X, uint8 (1)) %!test -%! X = dither (repmat(ones(3)/2,1,1,3), [0 0 0; 1 1 1], 4, 4); -%! assert (X, uint8([1 0 1; 0 1 0; 1 0 1])) +%! X = dither (repmat (ones (3)/2,1,1,3), [0 0 0; 1 1 1], 4, 4); +%! assert (X, uint8 ([1 0 1; 0 1 0; 1 0 1])) %!test -%! X = dither (repmat(ones(3)/4,1,1,3), [0 0 0; 1 1 1], 4, 4); -%! assert (X, uint8([0 0 0; 0 1 0; 0 0 0])) +%! X = dither (repmat (ones (3)/4,1,1,3), [0 0 0; 1 1 1], 4, 4); +%! assert (X, uint8 ([0 0 0; 0 1 0; 0 0 0])) %!test -%! X = dither (repmat(ones(3)*3/4,1,1,3), [0 0 0; 1 1 1], 4, 4); -%! assert (X, uint8([1 1 1; 1 0 1; 1 1 1])) +%! X = dither (repmat (ones (3)*3/4,1,1,3), [0 0 0; 1 1 1], 4, 4); +%! assert (X, uint8 ([1 1 1; 1 0 1; 1 1 1]))