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.

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:

1. Take the original MATLAB implementation from the paper
2. Run it in Octave with test images
3. Run your JavaScript implementation with the same images
4. 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]);
}

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.

The MATLAB Reference

The original implementation is 35 lines:

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:

1. Downsample both images by 2x using an averaging filter
2. Compute gradient magnitudes with Prewitt operators (3x3 edge detection)
3. 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²).

In MATLAB, this is a one-liner with conv2. In JavaScript, we loop through each pixel and its 8 neighbors:

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);
}

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]);
}

Breaking this down:

1. addpath — tells Octave where to find GMSD.m (the original MATLAB file from the paper)
2. imread — loads the test images
3. rgb2gray — converts to grayscale if needed (GMSD works on luminance only)
4. double() — MATLAB/Octave need floating-point input, not uint8
5. fprintf('%.15f') — outputs 15 decimal places, enough precision to catch subtle bugs

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 <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

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

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.

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.

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.

I ran the comparison. The results were surprising:

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);

The algorithm:

1. Downsample if the image is larger than 256px on any side
2. Apply a Gaussian window (11x11, σ=1.5) to compute local statistics
3. Calculate mean (μ), variance (σ²), and covariance for each window 4
4. Combine into the SSIM formula with stability constants C1 and C2
5. 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:

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:

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 < height; y++) {
    const rowStart = y * width;
    for (let x = pad; x < width - pad; x++) {
      let sum = 0;
      const srcStart = rowStart + x - pad;
      for (let k = 0; k < kernelSize; k++) {
        sum += input[srcStart + k] * kernel1d[k];
      }
      temp[rowStart + x] = sum;
    }
  }

  // Pass 2: Vertical convolution (temp -> 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 };
}

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]);
}

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);
});

Results

Compared to ssim.js on the same images:

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.

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]

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]

This matters when you translate indexing. MATLAB's img(row, col) becomes img[row * width + col] in JavaScript, not 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(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

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.

Color Space Coefficients

Converting RGB to grayscale seems simple. Multiply by coefficients, sum the channels. But which coefficients?

MATLAB's rgb2gray uses BT.601:

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

Others use the simple average:

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.

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.

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;
}

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.

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;
  }
}

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
  }
}

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<string, Float32Array>();

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 < height; y++) {
  for (let x = 0; x < width; x++) {
    const left = x > 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:

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.

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.

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.

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.

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.

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.