src/color/kmeans/index.js
function fillMatrixColumn(
  iMin,
  iMax,
  cluster,
  matrix,
  backtrackMatrix,
  sums,
  sumsOfSquares
) {
  if (iMin > iMax) {
    return;
  }

  // Start at midpoint between iMin and iMax
  let i = Math.floor((iMin + iMax) / 2);

  matrix[cluster][i] = matrix[cluster - 1][i - 1];
  backtrackMatrix[cluster][i] = i;

  let jlow = cluster; // the lower end for j

  if (iMin > cluster) {
    jlow = Math.max(jlow, backtrackMatrix[cluster][iMin - 1] || 0);
  }
  jlow = Math.max(jlow, backtrackMatrix[cluster - 1][i] || 0);

  let jhigh = i - 1; // the upper end for j
  if (iMax < matrix.length - 1) {
    jhigh = Math.min(jhigh, backtrackMatrix[cluster][iMax + 1] || 0);
  }

  let sji;
  let sjlowi;
  let ssqjlow;
  let ssqj;
  for (let j = jhigh; j >= jlow; --j) {
    sji = ssq(j, i, sums, sumsOfSquares);

    if (sji + matrix[cluster - 1][jlow - 1] >= matrix[cluster][i]) {
      break;
    }

    // Examine the lower bound of the cluster border
    sjlowi = ssq(jlow, i, sums, sumsOfSquares);

    ssqjlow = sjlowi + matrix[cluster - 1][jlow - 1];

    if (ssqjlow < matrix[cluster][i]) {
      // Shrink the lower bound
      matrix[cluster][i] = ssqjlow;
      backtrackMatrix[cluster][i] = jlow;
    }
    jlow++;

    ssqj = sji + matrix[cluster - 1][j - 1];
    if (ssqj < matrix[cluster][i]) {
      matrix[cluster][i] = ssqj;
      backtrackMatrix[cluster][i] = j;
    }
  }

  fillMatrixColumn(
    iMin,
    i - 1,
    cluster,
    matrix,
    backtrackMatrix,
    sums,
    sumsOfSquares
  );
  fillMatrixColumn(
    i + 1,
    iMax,
    cluster,
    matrix,
    backtrackMatrix,
    sums,
    sumsOfSquares
  );
}
function fillMatrixColumn(
  iMin,
  iMax,
  cluster,
  matrix,
  backtrackMatrix,
  sums,
  sumsOfSquares
) {
  if (iMin > iMax) {
    return;
  }

  // Start at midpoint between iMin and iMax
  let i = Math.floor((iMin + iMax) / 2);

  matrix[cluster][i] = matrix[cluster - 1][i - 1];
  backtrackMatrix[cluster][i] = i;

  let jlow = cluster; // the lower end for j

  if (iMin > cluster) {
    jlow = Math.max(jlow, backtrackMatrix[cluster][iMin - 1] || 0);
  }
  jlow = Math.max(jlow, backtrackMatrix[cluster - 1][i] || 0);

  let jhigh = i - 1; // the upper end for j
  if (iMax < matrix.length - 1) {
    jhigh = Math.min(jhigh, backtrackMatrix[cluster][iMax + 1] || 0);
  }

  let sji;
  let sjlowi;
  let ssqjlow;
  let ssqj;
  for (let j = jhigh; j >= jlow; --j) {
    sji = ssq(j, i, sums, sumsOfSquares);

    if (sji + matrix[cluster - 1][jlow - 1] >= matrix[cluster][i]) {
      break;
    }

    // Examine the lower bound of the cluster border
    sjlowi = ssq(jlow, i, sums, sumsOfSquares);

    ssqjlow = sjlowi + matrix[cluster - 1][jlow - 1];

    if (ssqjlow < matrix[cluster][i]) {
      // Shrink the lower bound
      matrix[cluster][i] = ssqjlow;
      backtrackMatrix[cluster][i] = jlow;
    }
    jlow++;

    ssqj = sji + matrix[cluster - 1][j - 1];
    if (ssqj < matrix[cluster][i]) {
      matrix[cluster][i] = ssqj;
      backtrackMatrix[cluster][i] = j;
    }
  }

  fillMatrixColumn(
    iMin,
    i - 1,
    cluster,
    matrix,
    backtrackMatrix,
    sums,
    sumsOfSquares
  );
  fillMatrixColumn(
    i + 1,
    iMax,
    cluster,
    matrix,
    backtrackMatrix,
    sums,
    sumsOfSquares
  );
}
function fillMatrixColumn ( iMin , iMax , cluster , matrix , backtrackMatrix , sums , sumsOfSquares ) { if ( iMin > iMax ) { return ; } // Start at midpoint between iMin and iMax let i = Math . floor ( ( iMin + iMax ) / 2 ) ; matrix [ cluster ] [ i ] = matrix [ cluster - 1 ] [ i - 1 ] ; backtrackMatrix [ cluster ] [ i ] = i ; let jlow = cluster ; // the lower end for j if ( iMin > cluster ) { jlow = Math . max ( jlow , backtrackMatrix [ cluster ] [ iMin - 1 ] || 0 ) ; } jlow = Math . max ( jlow , backtrackMatrix [ cluster - 1 ] [ i ] || 0 ) ; let jhigh = i - 1 ; // the upper end for j if ( iMax < matrix . length - 1 ) { jhigh = Math . min ( jhigh , backtrackMatrix [ cluster ] [ iMax + 1 ] || 0 ) ; } let sji ; let sjlowi ; let ssqjlow ; let ssqj ; for ( let j = jhigh ; j >= jlow ; -- j ) { sji = ssq ( j , i , sums , sumsOfSquares ) ; if ( sji + matrix [ cluster - 1 ] [ jlow - 1 ] >= matrix [ cluster ] [ i ] ) { break ; } // Examine the lower bound of the cluster border sjlowi = ssq ( jlow , i , sums , sumsOfSquares ) ; ssqjlow = sjlowi + matrix [ cluster - 1 ] [ jlow - 1 ] ; if ( ssqjlow < matrix [ cluster ] [ i ] ) { // Shrink the lower bound matrix [ cluster ] [ i ] = ssqjlow ; backtrackMatrix [ cluster ] [ i ] = jlow ; } jlow ++ ; ssqj = sji + matrix [ cluster - 1 ] [ j - 1 ] ; if ( ssqj < matrix [ cluster ] [ i ] ) { matrix [ cluster ] [ i ] = ssqj ; backtrackMatrix [ cluster ] [ i ] = j ; } } fillMatrixColumn ( iMin , i - 1 , cluster , matrix , backtrackMatrix , sums , sumsOfSquares ) ; fillMatrixColumn ( i + 1 , iMax , cluster , matrix , backtrackMatrix , sums , sumsOfSquares ) ; }
Function that recursively divides and conquers computations
for cluster j

@private
@param {number} iMin Minimum index in cluster to be computed
@param {number} iMax Maximum index in cluster to be computed
@param {number} cluster Index of the cluster currently being computed
@param {Array<Array<number>>} matrix
@param {Array<Array<number>>} backtrackMatrix
@param {Array<number>} sums
@param {Array<number>} sumsOfSquares
Function that recursively divides and conquers computations for cluster j