Porting Scientific Algorithms from MATLAB to JavaScript

Problem Most scientific algorithms live in MATLAB - image processing and analysis, machine learning fundamentals, etc. Researchers write MATLAB, publish papers, and share .m files. MATLAB .m But production code based on these algorithms runs elsewhere. If you're building a web app that needs image quality metrics, you're writing JavaScript. So you find a paper, read the algorithm, and reimplement it. Here's the problem: how do you know your implementation is correct? You can write unit tests. Check that identical images return 1.0, and that different images return a lower value. But that only catches obvious bugs. The subtle ones are wrong boundary handling, slightly off coefficients, or missing normalization steps. Your code runs and returns plausible numbers, but it's not doing what the paper describes. Existing JavaScript libraries have this problem. They're porting other low-level language ports, each generation drifting further from the original. Nobody validates against the MATLAB reference because MATLAB costs money, and running it in CI isn’t trivial. I ran into this when working on image comparison tools. I wanted reliable implementations of SSIM (Structural Similarity Index) and GMSD (Gradient Magnitude Similarity Deviation). The JavaScript options were either slow, inaccurate, or both. So I ported them myself and built a validation system to prove they match MATLAB. Idea: Use Octave for Verification GNU Octave is a free, open-source alternative to MATLAB. It runs the same .m files, produces the same outputs, and installs on any CI server. The idea is simple: GNU Octave .m Take the original MATLAB implementation from the paper Run it in Octave with test images Run your JavaScript implementation with the same images Compare the results Take the original MATLAB implementation from the paper Run it in Octave with test images Run your JavaScript implementation with the same images Compare the results This isn't unit testing. It's ground-truth validation. You're not checking that your code behaves correctly according to your understanding. You're checking that it matches the reference implementation exactly. Here's what that looks like in practice: function runMatlabGmsd(img1Path: string, img2Path: string): number { const matlabScript = ` addpath('${matlabPath}'); img1 = imread('${img1Path}'); img2 = imread('${img2Path}'); if size(img1, 3) == 3, img1 = rgb2gray(img1); end if size(img2, 3) == 3, img2 = rgb2gray(img2); end result = GMSD(double(img1), double(img2)); fprintf('%.15f', result); `; const output = execSync(`octave --eval "${matlabScript}"`, { encoding: 'utf-8', }); const match = output.match(/(\d+\.\d+)/); return parseFloat(match[0]); } function runMatlabGmsd(img1Path: string, img2Path: string): number { const matlabScript = ` addpath('${matlabPath}'); img1 = imread('${img1Path}'); img2 = imread('${img2Path}'); if size(img1, 3) == 3, img1 = rgb2gray(img1); end if size(img2, 3) == 3, img2 = rgb2gray(img2); end result = GMSD(double(img1), double(img2)); fprintf('%.15f', result); `; const output = execSync(`octave --eval "${matlabScript}"`, { encoding: 'utf-8', }); const match = output.match(/(\d+\.\d+)/); return parseFloat(match[0]); } You call Octave from Node.js, parse the floating-point result, and compare it to your JavaScript output. Not compromising ambiguity - taking 15 decimal places. Starting Simple: GMSD GMSD measures image similarity using gradients. It computes edge information for both images, compares them pixel-by-pixel, and returns a single score. Lower means more similar. Zero means identical. The algorithm is based on Xue, W., Zhang, L., Mou, X., & Bovik, A. C. (2013). "Gradient Magnitude Similarity Deviation: A Highly Efficient Perceptual Image Quality Index." IEEE Transactions on Image Processing. It's straightforward, well-documented, and has no existing JavaScript implementation. A clean slate to demonstrate the porting process. Xue, W., Zhang, L., Mou, X., & Bovik, A. C. (2013). "Gradient Magnitude Similarity Deviation: A Highly Efficient Perceptual Image Quality Index." IEEE Transactions on Image Processing The MATLAB Reference The original implementation is 35 lines: The original implementation function [score, quality_map] = GMSD(Y1, Y2) % Gradient Magnitude Similarity Deviation % Wufeng Xue, Lei Zhang, Xuanqin Mou, and Alan C. Bovik T = 170; Down_step = 2; dx = [1 0 -1; 1 0 -1; 1 0 -1]/3; dy = dx'; aveKernel = fspecial('average', 2); aveY1 = conv2(Y1, aveKernel, 'same'); aveY2 = conv2(Y2, aveKernel, 'same'); Y1 = aveY1(1:Down_step:end, 1:Down_step:end); Y2 = aveY2(1:Down_step:end, 1:Down_step:end); IxY1 = conv2(Y1, dx, 'same'); IyY1 = conv2(Y1, dy, 'same'); gradientMap1 = sqrt(IxY1.^2 + IyY1.^2); IxY2 = conv2(Y2, dx, 'same'); IyY2 = conv2(Y2, dy, 'same'); gradientMap2 = sqrt(IxY2.^2 + IyY2.^2); quality_map = (2*gradientMap1.*gradientMap2 + T) ./ (gradientMap1.^2 + gradientMap2.^2 + T); score = std2(quality_map); function [score, quality_map] = GMSD(Y1, Y2) % Gradient Magnitude Similarity Deviation % Wufeng Xue, Lei Zhang, Xuanqin Mou, and Alan C. Bovik T = 170; Down_step = 2; dx = [1 0 -1; 1 0 -1; 1 0 -1]/3; dy = dx'; aveKernel = fspecial('average', 2); aveY1 = conv2(Y1, aveKernel, 'same'); aveY2 = conv2(Y2, aveKernel, 'same'); Y1 = aveY1(1:Down_step:end, 1:Down_step:end); Y2 = aveY2(1:Down_step:end, 1:Down_step:end); IxY1 = conv2(Y1, dx, 'same'); IyY1 = conv2(Y1, dy, 'same'); gradientMap1 = sqrt(IxY1.^2 + IyY1.^2); IxY2 = conv2(Y2, dx, 'same'); IyY2 = conv2(Y2, dy, 'same'); gradientMap2 = sqrt(IxY2.^2 + IyY2.^2); quality_map = (2*gradientMap1.*gradientMap2 + T) ./ (gradientMap1.^2 + gradientMap2.^2 + T); score = std2(quality_map); The algorithm: The algorithm Downsample both images by 2x using an averaging filter Compute gradient magnitudes with Prewitt operators (3x3 edge detection) Calculate per-pixel similarity 4. Return the standard deviation of the similarity map Downsample both images by 2x using an averaging filter Compute gradient magnitudes with Prewitt operators (3x3 edge detection) Calculate per-pixel similarity 4. Return the standard deviation of the similarity map The JavaScript Port The core idea: images that look similar have similar edges. A photo and its slightly compressed version have edges in the same places. A photo and a completely different photo don't. GMSD captures this by computing "gradient magnitude" - essentially, how strong the edges are at each pixel. It uses the Prewitt operator, a classic 3x3 filter that detects horizontal and vertical intensity changes. The magnitude combines both directions: sqrt(horizontal² + vertical²). sqrt(horizontal² + vertical²) In MATLAB, this is a one-liner with conv2. In JavaScript, we loop through each pixel and its 8 neighbors: conv2 function computeGradientMagnitudes( luma: Float32Array, width: number, height: number ): Float32Array { const grad = new Float32Array(width * height); for (let y = 1; y < height - 1; y++) { for (let x = 1; x < width - 1; x++) { const idx = y * width + x; // Fetch 3x3 neighborhood const tl = luma[(y - 1) * width + (x - 1)]; const tc = luma[(y - 1) * width + x]; const tr = luma[(y - 1) * width + (x + 1)]; const ml = luma[y * width + (x - 1)]; const mr = luma[y * width + (x + 1)]; const bl = luma[(y + 1) * width + (x - 1)]; const bc = luma[(y + 1) * width + x]; const br = luma[(y + 1) * width + (x + 1)]; // Prewitt Gx = [1 0 -1; 1 0 -1; 1 0 -1]/3 const gx = (tl + ml + bl - tr - mr - br) / 3; // Prewitt Gy = [1 1 1; 0 0 0; -1 -1 -1]/3 const gy = (tl + tc + tr - bl - bc - br) / 3; grad[idx] = Math.sqrt(gx * gx + gy * gy); } } return grad; } function computeGradientMagnitudes( luma: Float32Array, width: number, height: number ): Float32Array { const grad = new Float32Array(width * height); for (let y = 1; y < height - 1; y++) { for (let x = 1; x < width - 1; x++) { const idx = y * width + x; // Fetch 3x3 neighborhood const tl = luma[(y - 1) * width + (x - 1)]; const tc = luma[(y - 1) * width + x]; const tr = luma[(y - 1) * width + (x + 1)]; const ml = luma[y * width + (x - 1)]; const mr = luma[y * width + (x + 1)]; const bl = luma[(y + 1) * width + (x - 1)]; const bc = luma[(y + 1) * width + x]; const br = luma[(y + 1) * width + (x + 1)]; // Prewitt Gx = [1 0 -1; 1 0 -1; 1 0 -1]/3 const gx = (tl + ml + bl - tr - mr - br) / 3; // Prewitt Gy = [1 1 1; 0 0 0; -1 -1 -1]/3 const gy = (tl + tc + tr - bl - bc - br) / 3; grad[idx] = Math.sqrt(gx * gx + gy * gy); } } return grad; } Once you have gradient magnitudes for both images, GMSD compares them pixel-by-pixel and returns the standard deviation of the differences. High deviation means the edges are in different places. Low deviation means a similar structure. function computeStdDev(values: Float32Array): number { const len = values.length; if (len === 0) return 0; let sum = 0; for (let i = 0; i < len; i++) { sum += values[i]; } const mean = sum / len; let variance = 0; for (let i = 0; i < len; i++) { const diff = values[i] - mean; variance += diff * diff; } variance /= len; return Math.sqrt(variance); } function computeStdDev(values: Float32Array): number { const len = values.length; if (len === 0) return 0; let sum = 0; for (let i = 0; i < len; i++) { sum += values[i]; } const mean = sum / len; let variance = 0; for (let i = 0; i < len; i++) { const diff = values[i] - mean; variance += diff * diff; } variance /= len; return Math.sqrt(variance); } Validation This is where Octave earns its keep. We run the original MATLAB code and our JavaScript implementation on the same images and compare the outputs. The test function wraps Octave in a shell call: function runMatlabGmsd(img1Path: string, img2Path: string): number { const matlabPath = join(__dirname, '../matlab'); const matlabScript = [ `addpath('${matlabPath}')`, `img1 = imread('${img1Path}')`, `img2 = imread('${img2Path}')`, `if size(img1, 3) == 3, img1 = rgb2gray(img1); end`, `if size(img2, 3) == 3, img2 = rgb2gray(img2); end`, `result = GMSD(double(img1), double(img2))`, `fprintf('%.15f', result)`, ].join('; '); const output = execSync(`octave --eval "${matlabScript}"`, { encoding: 'utf-8', }); const match = output.match(/(\d+\.\d+)/); return parseFloat(match[0]); } function runMatlabGmsd(img1Path: string, img2Path: string): number { const matlabPath = join(__dirname, '../matlab'); const matlabScript = [ `addpath('${matlabPath}')`, `img1 = imread('${img1Path}')`, `img2 = imread('${img2Path}')`, `if size(img1, 3) == 3, img1 = rgb2gray(img1); end`, `if size(img2, 3) == 3, img2 = rgb2gray(img2); end`, `result = GMSD(double(img1), double(img2))`, `fprintf('%.15f', result)`, ].join('; '); const output = execSync(`octave --eval "${matlabScript}"`, { encoding: 'utf-8', }); const match = output.match(/(\d+\.\d+)/); return parseFloat(match[0]); } Breaking this down: addpath — tells Octave where to find GMSD.m (the original MATLAB file from the paper) imread — loads the test images rgb2gray — converts to grayscale if needed (GMSD works on luminance only) double() — MATLAB/Octave need floating-point input, not uint8 fprintf('%.15f') — outputs 15 decimal places, enough precision to catch subtle bugs addpath — tells Octave where to find GMSD.m (the original MATLAB file from the paper) addpath GMSD.m imread — loads the test images imread rgb2gray — converts to grayscale if needed (GMSD works on luminance only) rgb2gray double() — MATLAB/Octave need floating-point input, not uint8 double() fprintf('%.15f') — outputs 15 decimal places, enough precision to catch subtle bugs fprintf('%.15f') The regex extracts the number from Octave's output. Octave prints some startup info, so we grab just the floating-point result. Now the actual test. Using Vitest, but Jest or any other framework works the same: import { execSync } from 'node:child_process'; import { describe, expect, it } from 'vitest'; import { gmsd } from './index'; describe('GMSD - MATLAB Comparison', () => { it('should match MATLAB reference', async () => { const tsResult = gmsd(img1.data, img2.data, undefined, width, height, { downsample: 1, c: 170, }); const matlabResult = runMatlabGmsd(img1Path, img2Path); const percentDiff = Math.abs(tsResult - matlabResult) / matlabResult * 100; console.log(`TypeScript: ${tsResult.toFixed(6)}`); console.log(`MATLAB: ${matlabResult.toFixed(6)}`); console.log(`Diff: ${percentDiff.toFixed(2)}%`); // Accept { it('should match MATLAB reference', async () => { const tsResult = gmsd(img1.data, img2.data, undefined, width, height, { downsample: 1, c: 170, }); const matlabResult = runMatlabGmsd(img1Path, img2Path); const percentDiff = Math.abs(tsResult - matlabResult) / matlabResult * 100; console.log(`TypeScript: ${tsResult.toFixed(6)}`); console.log(`MATLAB: ${matlabResult.toFixed(6)}`); console.log(`Diff: ${percentDiff.toFixed(2)}%`); // Accept <2% difference (boundary handling varies) expect(percentDiff).toBeLessThan(2); }); }); The key is setting the right tolerance. Too strict (0.0001%) and you'll chase floating-point ghosts. Too loose (10%) and real bugs slip through. I landed on 2% for GMSD after understanding the sources of the differences (more on that in Results). For CI, Octave installs in seconds on Ubuntu: name: Test on: [push, pull_request] jobs: test: runs-on: ubuntu-latest steps: - uses: actions/checkout@v4 - name: Install Octave run: sudo apt-get update && sudo apt-get install -y octave - name: Setup Node.js uses: actions/setup-node@v4 with: node-version: '20' - name: Install dependencies run: pnpm install - name: Run tests run: pnpm test name: Test on: [push, pull_request] jobs: test: runs-on: ubuntu-latest steps: - uses: actions/checkout@v4 - name: Install Octave run: sudo apt-get update && sudo apt-get install -y octave - name: Setup Node.js uses: actions/setup-node@v4 with: node-version: '20' - name: Install dependencies run: pnpm install - name: Run tests run: pnpm test Every push now validates your JavaScript against the MATLAB reference. If you accidentally break boundary handling or mess up a coefficient, the test fails. No more "it looks about right." Results Image Pair TypeScript MATLAB Difference 1a vs 1b 0.088934 0.089546 0.68% 2a vs 2b 0.142156 0.143921 1.23% 3a vs 3b 0.067823 0.068412 0.86% Image Pair TypeScript MATLAB Difference 1a vs 1b 0.088934 0.089546 0.68% 2a vs 2b 0.142156 0.143921 1.23% 3a vs 3b 0.067823 0.068412 0.86% Image Pair TypeScript MATLAB Difference Image Pair Image Pair TypeScript TypeScript MATLAB MATLAB Difference Difference 1a vs 1b 0.088934 0.089546 0.68% 1a vs 1b 1a vs 1b 0.088934 0.088934 0.089546 0.089546 0.68% 0.68% 2a vs 2b 0.142156 0.143921 1.23% 2a vs 2b 2a vs 2b 0.142156 0.142156 0.143921 0.143921 1.23% 1.23% 3a vs 3b 0.067823 0.068412 0.86% 3a vs 3b 3a vs 3b 0.067823 0.067823 0.068412 0.068412 0.86% 0.86% The 0.68-1.23% difference comes from boundary handling in conv2. MATLAB's 'same' mode handles image edges differently than my implementation. As a perceptual quality metric, this difference is acceptable because the relative rankings of image pairs remain identical. conv2 'same' Scaling Up: SSIM SSIM (Structural Similarity Index) is the industry standard for image quality assessment. It compares luminance, contrast, and structure between two images and returns a score from 0 to 1. Higher means more similar. 1.0 means identical. The algorithm comes from Wang, Z., Bovik, A. C., Sheikh, H. R., & Simoncelli, E. P (2004) IEEE Transactions on Image Processing paper. Image quality assessment: from error visibility to structural similarity. It's widely cited, well understood, and already has JavaScript implementations. The complication is we're not starting from scratch. We're entering a space with existing libraries, so we can directly compare accuracy. Wang, Z., Bovik, A. C., Sheikh, H. R., & Simoncelli, E. P (2004) IEEE Transactions on Image Processing The Existing Landscape The most popular JavaScript implementation is ssim.js. It works. It returns plausible numbers. But nobody has validated it against the MATLAB reference. ssim.js I ran the comparison. The results were surprising: Library Difference from MATLAB ssim.js 0.05%-0.73% @blazediff/ssim 0.00%-0.03% Library Difference from MATLAB ssim.js 0.05%-0.73% @blazediff/ssim 0.00%-0.03% Library Difference from MATLAB Library Library Difference from MATLAB Difference from MATLAB ssim.js 0.05%-0.73% ssim.js ssim.js 0.05%-0.73% 0.05%-0.73% @blazediff/ssim 0.00%-0.03% @blazediff/ssim @blazediff/ssim @blazediff/ssim 0.00%-0.03% 0.00%-0.03% That's up to 24x more accurate. Not because ssim.js is badly written — it's a reasonable implementation. But small decisions accumulate: slightly different Gaussian windows, missing auto-downsampling, different boundary handling. Each one adds an error. The MATLAB Reference Wang's reference implementation is more complex than GMSD. The key addition is automatic downsampling. Large images get scaled down before comparison. This matches human perception (we don't notice single-pixel differences in a 4K image) and improves performance. The core computation: % Automatic downsampling f = max(1, round(min(M,N)/256)); if(f > 1) lpf = ones(f,f); lpf = lpf/sum(lpf(:)); img1 = imfilter(img1, lpf, 'symmetric', 'same'); img2 = imfilter(img2, lpf, 'symmetric', 'same'); img1 = img1(1:f:end, 1:f:end); img2 = img2(1:f:end, 1:f:end); end C1 = (K(1)*L)^2; C2 = (K(2)*L)^2; window = window/sum(sum(window)); mu1 = filter2(window, img1, 'valid'); mu2 = filter2(window, img2, 'valid'); mu1_sq = mu1.*mu1; mu2_sq = mu2.*mu2; mu1_mu2 = mu1.*mu2; sigma1_sq = filter2(window, img1.*img1, 'valid') - mu1_sq; sigma2_sq = filter2(window, img2.*img2, 'valid') - mu2_sq; sigma12 = filter2(window, img1.*img2, 'valid') - mu1_mu2; ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2)) ./ ... ((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2)); mssim = mean2(ssim_map); % Automatic downsampling f = max(1, round(min(M,N)/256)); if(f > 1) lpf = ones(f,f); lpf = lpf/sum(lpf(:)); img1 = imfilter(img1, lpf, 'symmetric', 'same'); img2 = imfilter(img2, lpf, 'symmetric', 'same'); img1 = img1(1:f:end, 1:f:end); img2 = img2(1:f:end, 1:f:end); end C1 = (K(1)*L)^2; C2 = (K(2)*L)^2; window = window/sum(sum(window)); mu1 = filter2(window, img1, 'valid'); mu2 = filter2(window, img2, 'valid'); mu1_sq = mu1.*mu1; mu2_sq = mu2.*mu2; mu1_mu2 = mu1.*mu2; sigma1_sq = filter2(window, img1.*img1, 'valid') - mu1_sq; sigma2_sq = filter2(window, img2.*img2, 'valid') - mu2_sq; sigma12 = filter2(window, img1.*img2, 'valid') - mu1_mu2; ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2)) ./ ... ((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2)); mssim = mean2(ssim_map); The algorithm: The algorithm: Downsample if the image is larger than 256px on any side Apply a Gaussian window (11x11, σ=1.5) to compute local statistics Calculate mean (μ), variance (σ²), and covariance for each window 4 Combine into the SSIM formula with stability constants C1 and C2 Return the mean of all local SSIM values Downsample if the image is larger than 256px on any side Apply a Gaussian window (11x11, σ=1.5) to compute local statistics Calculate mean (μ), variance (σ²), and covariance for each window 4 Combine into the SSIM formula with stability constants C1 and C2 Return the mean of all local SSIM values The JavaScript Port Three parts needed careful attention: Gaussian windows, separable convolution, and the downsampling filter. Gaussian Windows MATLAB's fspecial('gaussian', 11, 1.5) creates an 11x11 Gaussian kernel. The JavaScript equivalent: fspecial('gaussian', 11, 1.5) function createGaussianWindow1D(size: number, sigma: number): Float32Array { const window = new Float32Array(size); const center = (size - 1) / 2; const twoSigmaSquared = 2 * sigma * sigma; let sum = 0; for (let i = 0; i < size; i++) { const d = i - center; const value = Math.exp(-(d * d) / twoSigmaSquared); window[i] = value; sum += value; } // Normalize so weights sum to 1 for (let i = 0; i < size; i++) { window[i] /= sum; } return window; } function createGaussianWindow1D(size: number, sigma: number): Float32Array { const window = new Float32Array(size); const center = (size - 1) / 2; const twoSigmaSquared = 2 * sigma * sigma; let sum = 0; for (let i = 0; i < size; i++) { const d = i - center; const value = Math.exp(-(d * d) / twoSigmaSquared); window[i] = value; sum += value; } // Normalize so weights sum to 1 for (let i = 0; i < size; i++) { window[i] /= sum; } return window; } Notice that this creates a 1D window, not a 2D one. That's intentional because Gaussian filters are separable, meaning a 2D convolution can be split into two 1D passes. Same result, half the operations. Separable Convolution MATLAB's filter2 applies a 2D kernel. For an 11x11 window on a 1000x1000 image, that's 121 million multiplications. Separable convolution cuts it to 22 million: filter2 function convolveSeparable( input: Float32Array, output: Float32Array, temp: Float32Array, width: number, height: number, kernel1d: Float32Array, kernelSize: number ): void { const pad = Math.floor(kernelSize / 2); // Pass 1: Horizontal convolution (input -> temp) for (let y = 0; y output) const outWidth = width - kernelSize + 1; const outHeight = height - kernelSize + 1; let outIdx = 0; for (let y = 0; y < outHeight; y++) { for (let x = 0; x < outWidth; x++) { let sum = 0; const srcX = x + pad; for (let k = 0; k < kernelSize; k++) { sum += temp[(y + k) * width + srcX] * kernel1d[k]; } output[outIdx++] = sum; } } } function convolveSeparable( input: Float32Array, output: Float32Array, temp: Float32Array, width: number, height: number, kernel1d: Float32Array, kernelSize: number ): void { const pad = Math.floor(kernelSize / 2); // Pass 1: Horizontal convolution (input -> temp) for (let y = 0; y output) const outWidth = width - kernelSize + 1; const outHeight = height - kernelSize + 1; let outIdx = 0; for (let y = 0; y < outHeight; y++) { for (let x = 0; x < outWidth; x++) { let sum = 0; const srcX = x + pad; for (let k = 0; k < kernelSize; k++) { sum += temp[(y + k) * width + srcX] * kernel1d[k]; } output[outIdx++] = sum; } } } Automatic Downsampling This is the part most JavaScript implementations skip. MATLAB's reference downsamples images larger than 256px using a box filter with symmetric padding: function downsampleImages( img1: Float32Array, img2: Float32Array, width: number, height: number, f: number ): { img1: Float32Array; img2: Float32Array; width: number; height: number } { // Create 1D averaging filter const filter1d = new Float32Array(f); const filterValue = 1 / f; for (let i = 0; i < f; i++) { filter1d[i] = filterValue; } // Apply separable filter with symmetric padding const temp = new Float32Array(width * height); const filtered1 = new Float32Array(width * height); const filtered2 = new Float32Array(width * height); convolveSeparableSymmetric(img1, filtered1, temp, width, height, filter1d, f); convolveSeparableSymmetric(img2, filtered2, temp, width, height, filter1d, f); // Subsample by taking every f-th pixel const newWidth = Math.floor(width / f); const newHeight = Math.floor(height / f); const downsampled1 = new Float32Array(newWidth * newHeight); const downsampled2 = new Float32Array(newWidth * newHeight); for (let y = 0; y < newHeight; y++) { for (let x = 0; x < newWidth; x++) { downsampled1[y * newWidth + x] = filtered1[y * f * width + x * f]; downsampled2[y * newWidth + x] = filtered2[y * f * width + x * f]; } } return { img1: downsampled1, img2: downsampled2, width: newWidth, height: newHeight }; } function downsampleImages( img1: Float32Array, img2: Float32Array, width: number, height: number, f: number ): { img1: Float32Array; img2: Float32Array; width: number; height: number } { // Create 1D averaging filter const filter1d = new Float32Array(f); const filterValue = 1 / f; for (let i = 0; i < f; i++) { filter1d[i] = filterValue; } // Apply separable filter with symmetric padding const temp = new Float32Array(width * height); const filtered1 = new Float32Array(width * height); const filtered2 = new Float32Array(width * height); convolveSeparableSymmetric(img1, filtered1, temp, width, height, filter1d, f); convolveSeparableSymmetric(img2, filtered2, temp, width, height, filter1d, f); // Subsample by taking every f-th pixel const newWidth = Math.floor(width / f); const newHeight = Math.floor(height / f); const downsampled1 = new Float32Array(newWidth * newHeight); const downsampled2 = new Float32Array(newWidth * newHeight); for (let y = 0; y < newHeight; y++) { for (let x = 0; x < newWidth; x++) { downsampled1[y * newWidth + x] = filtered1[y * f * width + x * f]; downsampled2[y * newWidth + x] = filtered2[y * f * width + x * f]; } } return { img1: downsampled1, img2: downsampled2, width: newWidth, height: newHeight }; } Skipping this step doesn't break anything obvious. You still get SSIM scores. But they won't match the reference, and the algorithm will run slower on large images. Validation Same pattern as GMSD. Call Octave, compare results: function runMatlabSsim(img1Path: string, img2Path: string): number { const matlabPath = join(__dirname, '../matlab'); const matlabScript = [ `addpath('${matlabPath}')`, `img1 = imread('${img1Path}')`, `img2 = imread('${img2Path}')`, `if size(img1, 3) == 3, img1 = rgb2gray(img1); end`, `if size(img2, 3) == 3, img2 = rgb2gray(img2); end`, `result = ssim(double(img1), double(img2))`, `fprintf('%.15f', result)`, ].join('; '); const output = execSync(`octave --eval "${matlabScript}"`, { encoding: 'utf-8', }); const match = output.match(/(\d+\.\d+)/); return parseFloat(match[0]); } function runMatlabSsim(img1Path: string, img2Path: string): number { const matlabPath = join(__dirname, '../matlab'); const matlabScript = [ `addpath('${matlabPath}')`, `img1 = imread('${img1Path}')`, `img2 = imread('${img2Path}')`, `if size(img1, 3) == 3, img1 = rgb2gray(img1); end`, `if size(img2, 3) == 3, img2 = rgb2gray(img2); end`, `result = ssim(double(img1), double(img2))`, `fprintf('%.15f', result)`, ].join('; '); const output = execSync(`octave --eval "${matlabScript}"`, { encoding: 'utf-8', }); const match = output.match(/(\d+\.\d+)/); return parseFloat(match[0]); } The tolerance is tighter here. SSIM is a well-defined algorithm with exact coefficients. If we're more than 0.05% off, something's wrong: it('should match MATLAB for images 1a vs 1b', async () => { const tsResult = ssim(png1.data, png2.data, undefined, width, height); const matlabResult = runMatlabSsim(img1Path, img2Path); const percentDiff = Math.abs(tsResult - matlabResult) / matlabResult * 100; // Strict: must be within 0.05% expect(percentDiff).toBeLessThan(0.05); }); it('should match MATLAB for images 1a vs 1b', async () => { const tsResult = ssim(png1.data, png2.data, undefined, width, height); const matlabResult = runMatlabSsim(img1Path, img2Path); const percentDiff = Math.abs(tsResult - matlabResult) / matlabResult * 100; // Strict: must be within 0.05% expect(percentDiff).toBeLessThan(0.05); }); Results Image Pair @blazediff/ssim MATLAB Difference 1a vs 1b 0.968753 0.968755 0.00% 2a vs 2b 0.912847 0.912872 0.00% 3a vs 3b 0.847621 0.847874 0.03% Image Pair @blazediff/ssim MATLAB Difference 1a vs 1b 0.968753 0.968755 0.00% 2a vs 2b 0.912847 0.912872 0.00% 3a vs 3b 0.847621 0.847874 0.03% Image Pair @blazediff/ssim MATLAB Difference Image Pair Image Pair @blazediff/ssim @blazediff/ssim MATLAB MATLAB Difference Difference 1a vs 1b 0.968753 0.968755 0.00% 1a vs 1b 1a vs 1b 0.968753 0.968753 0.968755 0.968755 0.00% 0.00% 2a vs 2b 0.912847 0.912872 0.00% 2a vs 2b 2a vs 2b 0.912847 0.912847 0.912872 0.912872 0.00% 0.00% 3a vs 3b 0.847621 0.847874 0.03% 3a vs 3b 3a vs 3b 0.847621 0.847621 0.847874 0.847874 0.03% 0.03% Compared to ssim.js on the same images: Image Pair ssim.js MATLAB Difference 1a vs 1b 0.969241 0.968755 0.05% 2a vs 2b 0.916523 0.912872 0.40% 3a vs 3b 0.853842 0.847874 0.73% Image Pair ssim.js MATLAB Difference 1a vs 1b 0.969241 0.968755 0.05% 2a vs 2b 0.916523 0.912872 0.40% 3a vs 3b 0.853842 0.847874 0.73% Image Pair ssim.js MATLAB Difference Image Pair Image Pair ssim.js ssim.js MATLAB MATLAB Difference Difference 1a vs 1b 0.969241 0.968755 0.05% 1a vs 1b 1a vs 1b 0.969241 0.969241 0.968755 0.968755 0.05% 0.05% 2a vs 2b 0.916523 0.912872 0.40% 2a vs 2b 2a vs 2b 0.916523 0.916523 0.912872 0.912872 0.40% 0.40% 3a vs 3b 0.853842 0.847874 0.73% 3a vs 3b 3a vs 3b 0.853842 0.853842 0.847874 0.847874 0.73% 0.73% The 0.00-0.03% difference in my implementation is due to floating-point rounding. It’s unavoidable when translating between languages. The scores are functionally identical to MATLAB. The ssim.js differences come from algorithmic choices: no auto-downsampling, slightly different Gaussian approximation, and different boundary handling. None are bugs, but they compound into measurable error. Common Pitfalls Porting MATLAB to JavaScript looks straightforward until it isn't. These are the issues that cost me time. Array Indexing MATLAB arrays start at 1. JavaScript arrays start at 0. Everyone knows this. Everyone still gets bitten by it. The bug is rarely obvious. You won't get an off-by-one error on the first pixel. You'll get slightly wrong values at image boundaries, where loops like for i = 1:height translate to for (let i = 0; i < height; i++) but the neighbor access img(i-1, j) becomes img[(i-1) * width + j], which hits index -1 when i is 0. for i = 1:height for (let i = 0; i < height; i++) img(i-1, j) img[(i-1) * width + j] i Fix: draw out a 3x3 example on paper before writing the loop. Check your boundary conditions explicitly. Column-Major vs Row-Major MATLAB stores matrices column-by-column. JavaScript TypedArrays store row-by-row. A 3x3 matrix in MATLAB: [1 4 7] [2 5 8] → stored as [1, 2, 3, 4, 5, 6, 7, 8, 9] [3 6 9] [1 4 7] [2 5 8] → stored as [1, 2, 3, 4, 5, 6, 7, 8, 9] [3 6 9] The same matrix in JavaScript: [1 4 7] [2 5 8] → stored as [1, 4, 7, 2, 5, 8, 3, 6, 9] [3 6 9] [1 4 7] [2 5 8] → stored as [1, 4, 7, 2, 5, 8, 3, 6, 9] [3 6 9] This matters when you translate indexing. MATLAB's img(row, col) becomes img[row * width + col] in JavaScript, not img[col * height + row]. img(row, col) img[row * width + col] img[col * height + row] Most image libraries already hand you row-major data, so you're fine. But if you're copying MATLAB matrix operations verbatim, watch out. Boundary Handling MATLAB's conv2, filter2, and imfilter have different default behaviors: conv2 filter2 imfilter conv2(A, K) — no padding, output shrinks conv2(A, K, 'same') — zero padding, output same size as input filter2(K, A, 'valid') — no padding, output shrinks imfilter(A, K, 'symmetric') — mirror padding at edges conv2(A, K) — no padding, output shrinks conv2(A, K) conv2(A, K, 'same') — zero padding, output same size as input conv2(A, K, 'same') filter2(K, A, 'valid') — no padding, output shrinks filter2(K, A, 'valid') imfilter(A, K, 'symmetric') — mirror padding at edges imfilter(A, K, 'symmetric') Get this wrong, and your results will differ at every pixel near the border. For a 1000x1000 image with an 11x11 kernel, that's ~4% of your pixels. The SSIM reference uses imfilter with 'symmetric' padding for downsampling, then filter2 with 'valid' mode for the main computation. Miss either detail, and you'll wonder why your numbers are 2% off. imfilter 'symmetric' filter2 'valid' Color Space Coefficients Converting RGB to grayscale seems simple. Multiply by coefficients, sum the channels. But which coefficients? MATLAB's rgb2gray uses BT.601: rgb2gray Y = 0.298936 * R + 0.587043 * G + 0.114021 * B Y = 0.298936 * R + 0.587043 * G + 0.114021 * B Some JavaScript libraries use BT.709: Y = 0.2126 * R + 0.7152 * G + 0.0722 * B Y = 0.2126 * R + 0.7152 * G + 0.0722 * B Others use the simple average: Y = (R + G + B) / 3 Y = (R + G + B) / 3 The difference is negligible for most images. But if you're validating against MATLAB, use the exact BT.601 coefficients. Otherwise, you'll chase phantom bugs that are really just a grayscale conversion mismatch. MATLAB’s Implicit Behaviours MATLAB does things automatically that JavaScript won't: Auto-casting: MATLAB's double(img) converts uint8 (0-255) to float (0.0-255.0). If your JavaScript reads PNG data as Uint8ClampedArray and you forget to convert, your variance calculations will overflow. Auto-casting double(img) Default parameters: MATLAB functions have default values buried in their documentation. SSIM uses K = [0.01, 0.03] , L = 255, window size 11, sigma 1.5. Miss one, and your implementation diverges. Default parameters K = [0.01, 0.03] L = 255 Optimizing Your JavaScript Accuracy comes first. But once your implementation matches MATLAB, performance matters. Scientific algorithms process millions of pixels. Small inefficiencies multiply. These optimizations made my SSIM implementation 25-70% faster than ssim.js while maintaining accuracy. Use TypedArrays Regular JavaScript arrays are flexible. They can hold mixed types, grow dynamically, and have convenient methods. That flexibility costs performance. // Slow: regular array const pixels = []; for (let i = 0; i < width * height; i++) { pixels.push(image[i] * 0.299); } // Fast: TypedArray const pixels = new Float32Array(width * height); for (let i = 0; i < width * height; i++) { pixels[i] = image[i] * 0.299; } // Slow: regular array const pixels = []; for (let i = 0; i < width * height; i++) { pixels.push(image[i] * 0.299); } // Fast: TypedArray const pixels = new Float32Array(width * height); for (let i = 0; i < width * height; i++) { pixels[i] = image[i] * 0.299; } TypedArrays have fixed size, fixed type, and contiguous memory. The JavaScript engine knows precisely how to optimize them. For numerical code, this is a 2-5x speedup. Use Float32Array for most computations. Use Float64Array when you need extra precision (accumulating sums over large images). Use Uint8ClampedArray for final image output. Float32Array Float64Array Uint8ClampedArray Separable Filters A 2D convolution with an NxN kernel requires N² multiplications per pixel. For SSIM's 11x11 Gaussian, that's 121 operations per pixel. But Gaussian filters are separable. A 2D Gaussian is the outer product of two 1D Gaussians. Instead of one 11x11 pass, you do two 11x1 passes: 22 operations per pixel. // Slow: 2D convolution (N² operations per pixel) for (let y = 0; y < height; y++) { for (let x = 0; x < width; x++) { let sum = 0; for (let ky = 0; ky < kernelSize; ky++) { for (let kx = 0; kx < kernelSize; kx++) { sum += getPixel(x + kx, y + ky) * kernel2d[ky][kx]; } } output[y * width + x] = sum; } } // Fast: separable convolution (2N operations per pixel) // Pass 1: horizontal for (let y = 0; y < height; y++) { for (let x = 0; x < width; x++) { let sum = 0; for (let k = 0; k < kernelSize; k++) { sum += getPixel(x + k, y) * kernel1d[k]; } temp[y * width + x] = sum; } } // Pass 2: vertical for (let y = 0; y < height; y++) { for (let x = 0; x < width; x++) { let sum = 0; for (let k = 0; k < kernelSize; k++) { sum += temp[(y + k) * width + x] * kernel1d[k]; } output[y * width + x] = sum; } } // Slow: 2D convolution (N² operations per pixel) for (let y = 0; y < height; y++) { for (let x = 0; x < width; x++) { let sum = 0; for (let ky = 0; ky < kernelSize; ky++) { for (let kx = 0; kx < kernelSize; kx++) { sum += getPixel(x + kx, y + ky) * kernel2d[ky][kx]; } } output[y * width + x] = sum; } } // Fast: separable convolution (2N operations per pixel) // Pass 1: horizontal for (let y = 0; y < height; y++) { for (let x = 0; x < width; x++) { let sum = 0; for (let k = 0; k < kernelSize; k++) { sum += getPixel(x + k, y) * kernel1d[k]; } temp[y * width + x] = sum; } } // Pass 2: vertical for (let y = 0; y < height; y++) { for (let x = 0; x < width; x++) { let sum = 0; for (let k = 0; k < kernelSize; k++) { sum += temp[(y + k) * width + x] * kernel1d[k]; } output[y * width + x] = sum; } } This works for any separable kernel: Gaussian, box filter, or Sobel. Check if your kernel is separable before implementing 2D convolution. Reuse Buffers Memory allocation is expensive. Garbage collection is worse. Allocate once, reuse everywhere. // Slow: allocate per operation function computeVariance(img: Float32Array): Float32Array { const squared = new Float32Array(img.length); // allocation const filtered = new Float32Array(img.length); // allocation // ... } // Fast: pre-allocate and reuse class SSIMComputer { private temp: Float32Array; private squared: Float32Array; constructor(maxSize: number) { this.temp = new Float32Array(maxSize); this.squared = new Float32Array(maxSize); } compute(img: Float32Array): number { // reuse this.temp and this.squared } } // Slow: allocate per operation function computeVariance(img: Float32Array): Float32Array { const squared = new Float32Array(img.length); // allocation const filtered = new Float32Array(img.length); // allocation // ... } // Fast: pre-allocate and reuse class SSIMComputer { private temp: Float32Array; private squared: Float32Array; constructor(maxSize: number) { this.temp = new Float32Array(maxSize); this.squared = new Float32Array(maxSize); } compute(img: Float32Array): number { // reuse this.temp and this.squared } } For SSIM, I allocate all buffers upfront: grayscale images, squared images, filtered outputs, and the SSIM map. One allocation at the start, zero during computation. Cache Computed Values Some values get reused. Don’t compute them twice. // Slow: recompute Gaussian window every call function ssim(img1, img2, width, height) { const window = createGaussianWindow(11, 1.5); // expensive // ... } // Fast: cache by parameters const windowCache = new Map (); function getGaussianWindow(size: number, sigma: number): Float32Array { const key = `${size}_${sigma}`; let window = windowCache.get(key); if (!window) { window = createGaussianWindow(size, sigma); windowCache.set(key, window); } return window; } // Slow: recompute Gaussian window every call function ssim(img1, img2, width, height) { const window = createGaussianWindow(11, 1.5); // expensive // ... } // Fast: cache by parameters const windowCache = new Map (); function getGaussianWindow(size: number, sigma: number): Float32Array { const key = `${size}_${sigma}`; let window = windowCache.get(key); if (!window) { window = createGaussianWindow(size, sigma); windowCache.set(key, window); } return window; } The Gaussian window for SSIM is always 11x11 with σ=1.5. Computing it takes microseconds. But when you're processing thousands of images, microseconds add up. Avoid Bounds Checking in Hot Loops JavaScript arrays check bounds on every access. For interior pixels where you know indices are valid, this is wasted work. // Slow: bounds checking on every access for (let y = 0; y 0 ? img[y * width + x - 1] : 0; const right = x < width - 1 ? img[y * width + x + 1] : 0; // ... } } // Fast: handle borders separately // Interior pixels (no bounds checking needed) for (let y = 1; y 0 ? img[y * width + x - 1] : 0; const right = x < width - 1 ? img[y * width + x + 1] : 0; // ... } } // Fast: handle borders separately // Interior pixels (no bounds checking needed) for (let y = 1; y < height - 1; y++) { for (let x = 1; x < width - 1; x++) { const left = img[y * width + x - 1]; // always valid const right = img[y * width + x + 1]; // always valid // ... } } // Handle border pixels separately with bounds checking This is uglier code. But for a 1000x1000 image, you're removing millions of conditional checks. Results Benchmarked on the same images, same machine: Implementation Time (avg) vs MATLAB accuracy ssim.js 86ms 0.05-0.73% @blazediff/ssim 64ms 0.00-0.03% Implementation Time (avg) vs MATLAB accuracy ssim.js 86ms 0.05-0.73% @blazediff/ssim 64ms 0.00-0.03% Implementation Time (avg) vs MATLAB accuracy Implementation Implementation Time (avg) Time (avg) vs MATLAB accuracy vs MATLAB accuracy ssim.js 86ms 0.05-0.73% ssim.js ssim.js 86ms 86ms 0.05-0.73% 0.05-0.73% @blazediff/ssim 64ms 0.00-0.03% @blazediff/ssim @blazediff/ssim 64ms 64ms 0.00-0.03% 0.00-0.03% 25% faster and more accurate. Not because of clever algorithms — because of careful engineering. TypedArrays, separable filters, buffer reuse, cached windows. For the Hitchhiker's SSIM variant (uses integral images instead of convolution), the gap is wider: 70% faster than ssim.js on large images. Hitchhiker's SSIM Takeaways Porting scientific algorithms from MATLAB to JavaScript is mechanical once you have the proper process. Use the reference implementation. Papers describe algorithms. MATLAB code implements them. These aren't the same thing. Edge cases, default parameters, boundary handling - they're in the code, not the paper. Start with the .m file. Use the reference implementation .m Validate with Octave. It's free, runs everywhere, and executes MATLAB code without modification. Set up ground-truth tests early. Run them in CI. When your numbers match to 4 decimal places, you're done. When they don't, you know immediately. Validate with Octave Accuracy first, then performance. Get your implementation matching MATLAB before optimizing. A fast wrong answer is still wrong. Once you're accurate, the optimizations are straightforward: TypedArrays, separable filters, buffer reuse. Accuracy first, then performance Document your differences. Sometimes you won't match exactly. Boundary handling, floating-point precision, and intentional simplifications. That's fine. Know why you differ and by how much. Write it down. Document your differences The gap between academic MATLAB and production JavaScript is smaller than it looks. The algorithms are the same. The math is the same. You need a way to verify you got it right. The SSIM and GMSD implementations from this article are available at [blazediff.dev](https://blazediff.dev). MIT licensed, zero dependencies, validated against MATLAB.