This Image Processing Reference section displays the code for an example program that demonstrates how to calculate the singular value decomposition for the color channels of an image.

SvdChannels.php

<?php
  function Hypotenuse($dA, $dB) {
    $dAbsA = abs($dA);
    $dAbsB = abs($dB);
    if ($dAbsA > $dAbsB) {
      $dQ = ($dAbsB/$dAbsA);
      return $dAbsA*sqrt(1 + $dQ*$dQ);
    } else if ($dAbsB == 0.0) {
      return 0.0;
    } else {
      $dQ = ($dAbsA/$dAbsB);
      return $dAbsB*sqrt(1 + $dQ*$dQ);
    }
  }

  function Sign($dA, $dB) {
    return ($dB >= 0.0
      ? ($dA >= 0.0 ? $dA : -$dA)
      : ($dA >= 0.0 ? -$dA : $dA));
  }

  // A is the m by n matrix that is to be decomposed. The matrix A = UxWx(TransposeOfV)
  // U is returned in the matrix A. A and Vt are 1D arrays that are large enough to hold matrices.
  function SVD($daAtoU, $iM, $iN, $daW, $daVt) {
    // Allocate an array for calculations
    $daRV1 = new SplFixedArray($iN);
    $dNorm = 0.0;
    // These three are probably not needed outside the loop
    $dScale = 0.0;
    $dG = 0.0;
    $dS = 0.0;
    $iL = 0;
    $iNM = 0;
    $dZ = 0.0;
    // Householder reflections to reduce the matrix to a bidiagonal matrix.
    for ($i = 0; $i < $iN; ++$i) {
      // It looks like $iL is always $i + 1
      $iL         = $i + 1;
      $daRV1[$i]  = ($dG*$dScale);
      $dScale     = 0.0;
      $dG         = 0.0;
      $dS         = 0.0;
      if ($i < $iM) {
        for ($k = $i; $k < $iM; ++$k) {
          $dScale += abs($daAtoU[$k*$iN + $i]);
        }
        if ($dScale != 0.0) {
          for ($k = $i; $k < $iM; ++$k) {
            $daAtoU[$k*$iN + $i] /= $dScale;
            $dS += ($daAtoU[$k*$iN + $i]*$daAtoU[$k*$iN + $i]);
          }
          $dF = $daAtoU[$i*$iN + $i];
          // This seems to always be the opposite value. Check it! Perhaps, just $dG = -$dF
          $dG = -Sign(sqrt($dS), $dF);
          $dH = ($dF*$dG - $dS);
          $daAtoU[$i*$iN + $i] = ($dF - $dG);
          for ($j = $iL; $j < $iN; ++$j) {
            $dS = 0.0;
            for ($k = $i; $k < $iM; ++$k) {
              $dS += ($daAtoU[$k*$iN + $i]*$daAtoU[$k*$iN + $j]);
            }
            $dF = ($dS/$dH);
            for ($k = $i; $k < $iM; ++$k) {
              $daAtoU[$k*$iN + $j] += ($dF*$daAtoU[$k*$iN + $i]);
            }
          }
          for ($k = $i; $k < $iM; ++$k) {
            $daAtoU[$k*$iN + $i] *= $dScale;
          }
        }
      } // if ($i < $iM)
      $daW[$i] = ($dScale*$dG);
      $dScale  = 0.0;
      $dG      = 0.0;
      $dS      = 0.0;
      if (($i < $iM) && ($i != ($iN - 1))) {
        for ($k = $iL; $k < $iN; ++$k) {
          $dScale += abs($daAtoU[$i*$iN + $k]);
        }
        if ($dScale != 0.0) {
          for ($k = $iL; $k < $iN; ++$k) {
            $daAtoU[$i*$iN + $k] /= $dScale;
            $dS += ($daAtoU[$i*$iN + $k]*$daAtoU[$i*$iN + $k]);
          }
          $dF = $daAtoU[$i*$iN + $iL];
          // This seems to always be the opposite value. Check it! Perhaps, just $dG = -$dF
          $dG = -Sign(sqrt($dS), $dF);
          $dH = ($dF*$dG - $dS);
          $daAtoU[$i*$iN + $iL] = ($dF - $dG);
          for ($k = $iL; $k < $iN; ++$k) {
            $daRV1[$k] = ($daAtoU[$i*$iN + $k]/$dH);
          }
          for ($j = $iL; $j < $iM; ++$j) {
            $dS = 0.0;
            for ($k = $iL; $k < $iN; ++$k) {
              $dS += ($daAtoU[$j*$iN + $k]*$daAtoU[$i*$iN + $k]);
            }
            for ($k = $iL; $k < $iN; ++$k) {
              $daAtoU[$j*$iN + $k] += ($dS*$daRV1[$k]);
            }
          }
          for ($k = $iL; $k < $iN; ++$k) {
            $daAtoU[$i*$iN + $k] *= $dScale;
          }
        }
      } //if (($i < $iM) && ($i != ($iN - 1))) {
      $dNorm = max($dNorm, abs($daW[$i]) + abs($daRV1[$i]));
    }
    // Accumulation of right-hand transformations
    for ($i = $iN - 1; $i >= 0; --$i) {
      if ($i < ($iN - 1)) {
        if ($dG != 0.0) {
          for ($j = $iL; $j < $iN; ++$j) {
            $daVt[$j*$iN + $i] = (($daAtoU[$i*$iN + $j]/$daAtoU[$i*$iN + $iL])/$dG);
          }
          for ($j = $iL; $j < $iN; ++$j) {
            $dS = 0.0;
            for ($k = $iL; $k < $iN; ++$k) {
              $dS += ($daAtoU[$i*$iN + $k]*$daVt[$k*$iN + $j]);
            }
            for ($k = $iL; $k < $iN; ++$k) {
              $daVt[$k*$iN + $j] += ($dS*$daVt[$k*$iN + $i]);
            }
          }
        } // if ($dG > 0.0) {
        for ($j = $iL; $j < $iN; ++$j) {
          $daVt[$i*$iN + $j] = 0.0;
          $daVt[$j*$iN + $i] = 0.0;
        }
      }
      $daVt[$i*$iN + $i] = 1.0;
      $dG = $daRV1[$i];
      $iL = $i;
    }
    // Accumulation of left-hand transformations
    for ($i = IntVal(Min($iN-1, $iM-1)); $i >= 0; --$i) {
      $iL = $i + 1;
      $dG = $daW[$i];
      for ($j = $iL; $j < $iN; ++$j) {
        $daAtoU[$i*$iN + $j] = 0.0;
      }
      if ($dG != 0.0) {
        $dG = (1.0/$dG);
        for ($j = $iL; $j < $iN; ++$j) {
          $dS = 0.0;
          for ($k = $iL; $k < $iM; ++$k) {
            $dS += ($daAtoU[$k*$iN + $i]*$daAtoU[$k*$iN + $j]);
          }
          $dF = ($dS/$daAtoU[$i*$iN + $i])*$dG;
          for ($k = $i; $k < $iM; ++$k) {
            $daAtoU[$k*$iN + $j] += ($dF*$daAtoU[$k*$iN + $i]);
          }
        }
        for ($j = $i; $j < $iM; ++$j) {
          $daAtoU[$j*$iN + $i] *= $dG;
        }
      } else {
        for ($j = $i; $j < $iM; ++$j) {
          $daAtoU[$j*$iN + $i] = 0.0;
        }
      }
      $daAtoU[$i*$iN + $i] += 1.0;
    }
    // Diagonalization of the bidiagonal form
    for ($k = $iN - 1; $k >= 0; --$k) {
      for ($iIter = 0; $iIter < 30; ++$iIter) {
        $bFlag = true;
        // Test for splitting. Note: $daRV1[$i] is always zero
        // This is the while loop equivalent of the for loop in the code.
        $iL = $k;
        $bLoop = ($iL >= 0);
        // Test for splitting. $daRV1[0] is always zero
        while ($bLoop) {
          $iNM = $iL - 1;
          if ((abs($daRV1[$iL]) + $dNorm) == $dNorm) {
            $bFlag = FALSE;
            $bLoop = FALSE;
          } else if ((abs($daW[$iNM]) + $dNorm) == $dNorm) {
            $bLoop = FALSE;
          } else {
            --$iL;
            if ($iL < 0) {
              $bLoop = FALSE;
            }
          }
        }
        if ($bFlag) {
          // Cancellation of $daRV1[$iL], if L > 0
          $dC = 0.0;
          $dS = 1.0;
          for ($i = $iL; $i <= $k; ++$i) {
            $dF = ($dS*$daRV1[$i]);
            $daRV1[$i] = ($dC*$daRV1[$i]);
            if ((abs($dF)+$dNorm) == $dNorm) {
              break;
            }
            $dG = $daW[$i];
            $dH = Hypotenuse($dF, $dG);
            $daW[$i] = $dH;
            $dH = (1.0/$dH);
            $dC = ($dG*$dH);
            $dS = (-$dF*$dH);
            // Perform rotations
            for ($j = 0; $j < $iM; ++$j) {
              $dY = $daAtoU[$j*$iN + $iNM];
              $dZ = $daAtoU[$j*$iN + $i];
              $daAtoU[$j*$iN + $iNM] = (($dY*$dC) + ($dZ*$dS));
              $daAtoU[$j*$iN + $i] = (($dZ*$dC) - ($dY*$dS));
            }
          }
        } // if ($bFlag)
        $dZ = $daW[$k];
        // Convergence
        if ($iL == $k) {
          // Singular value is made nonnegative
          if ($dZ < 0.0) {
            $daW[$k] = -$dZ;
            for ($j = 0; $j < $iN; ++$j) {
              $daVt[$j*$iN + $k] = -$daVt[$j*$iN + $k];
            }
          }
          // Break out of the outer loop
          break;
        }
        if ($iIter == 29) { // No convergence. What should we do??
        }
        $dX = $daW[$iL];
        $iNM = $k - 1;
        $dY = $daW[$iNM];
        $dG = $daRV1[$iNM];
        $dH = $daRV1[$k];
        $dF = ((($dY-$dZ)*($dY+$dZ)+($dG-$dH)*($dG+$dH))/(2.0*$dH*$dY));
        $dG = Hypotenuse($dF, 1.0);
        $dF = ( ($dX-$dZ)*($dX+$dZ)+$dH*( ($dY/($dF+Sign($dG,$dF))) -$dH) )/$dX;
        $dC = 1.0;
        $dS = 1.0;
        // Next QR Transformation
        for ($j = $iL; $j <= $iNM; ++$j) {
          $i = $j + 1;
          $dG = $daRV1[$i];
          $dY = $daW[$i];
          $dH = $dS*$dG;
          $dG = $dC*$dG;
          $dZ = Hypotenuse($dF, $dH);
          $daRV1[$j] = $dZ;
          $dC = $dF/$dZ;
          $dS = $dH/$dZ;
          $dF = ($dX*$dC + $dG*$dS);
          $dG = ($dG*$dC - $dX*$dS);
          $dH = $dY*$dS;
          $dY *= $dC;
          for ($jj = 0; $jj < $iN; ++$jj) {
            $dX = $daVt[$jj*$iN + $j];
            $dZ = $daVt[$jj*$iN + $i];
            $daVt[$jj*$iN + $j] = $dX*$dC + $dZ*$dS;
            $daVt[$jj*$iN + $i] = $dZ*$dC - $dX*$dS;
          }
          $dZ = Hypotenuse($dF, $dH);
          $daW[$j] = $dZ;
          if ($dZ != 0.0) {
            $dZ = (1.0/$dZ);
            $dC = $dF*$dZ;
            $dS = $dH*$dZ;
          }
          $dF = ($dG*$dC + $dY*$dS);
          $dX = ($dY*$dC - $dG*$dS);
          for ($jj = 0; $jj < $iM; ++$jj) {
            $dY = $daAtoU[$jj*$iN + $j];
            $dZ = $daAtoU[$jj*$iN + $i];
            $daAtoU[$jj*$iN + $j] = ($dY*$dC + $dZ*$dS);
            $daAtoU[$jj*$iN + $i] = ($dZ*$dC - $dY*$dS);
          }
        }
        $daRV1[$iL] = 0.0;
        $daRV1[$k] = $dF;
        $daW[$k] = $dX;
      } //for ($iIter = 0; $iIter < 30; ++$iIter) {
    }
  } // End of function

  function Clamp($iGray) {
    if ($iGray < 0) {
      return 0;
    } elseif ($iGray > 255) {
      return 255;
    } else {
      return $iGray;
    }
  }

  function CreateSvdColorChannelsImage($qColorImage) {
    $iW   = ImageSX($qColorImage);
    $iH   = ImageSY($qColorImage);
    // Allocate the color channel images
    $qImageR   = ImageCreate($iW, 5*$iH);
    $qImageG   = ImageCreate($iW, 5*$iH);
    $qImageB   = ImageCreate($iW, 5*$iH);
    // This will store the full color combination image
    $qImageColor = ImageCreateTrueColor($iW, 5*$iH);
    for ($iC = 0; $iC < 256; ++$iC) {
      ImageColorAllocate($qImageR, $iC,   0,   0);
      ImageColorAllocate($qImageG,   0, $iC,   0);
      ImageColorAllocate($qImageB,   0,   0, $iC);
    }
    // Create the Red, Green, and Blue Channel images
    for ($j = 0; $j < $iH; ++$j) {
      for ($i = 0; $i < $iW; ++$i) {
        $qColor  = ImageColorAt($qColorImage, $i, $j);
        $iRed    = (($qColor >> 16) & 0xFF);
        $iGreen  = (($qColor >> 8) & 0xFF);
        $iBlue   = ($qColor & 0xFF);
        ImageSetPixel($qImageR, $i, $j, $iRed);
        ImageSetPixel($qImageG, $i, $j, $iGreen);
        ImageSetPixel($qImageB, $i, $j, $iBlue);
      }
    }
    $qaChannelImages = array();
    $qaChannelImages[0] = $qImageR;
    $qaChannelImages[1] = $qImageG;
    $qaChannelImages[2] = $qImageB;

    // We will reuse these arrays for each channel
    $daA	 = new SplFixedArray($iW*$iH);
    $daW	 = new SplFixedArray($iW);
    $daVt	 = new SplFixedArray($iW*$iW);
    $daAccum = new SplFixedArray($iW*$iH);

    for ($iChannel = 0; $iChannel < 3; ++$iChannel) {
      // Copy the color image to the array A
      for ($j = 0; $j < $iH; ++$j) {
        for ($i = 0; $i < $iW; ++$i) {
          // This is a number  between 0 and 255
          $iPixel  = ImageColorAt($qaChannelImages[$iChannel], $i, $j);
          $daA[$j*$iW + $i] = $iPixel;
          $daAccum[$j*$iW + $i] = 0.0;
        }
      }
      // Call the SVD routine
      SVD($daA, $iH, $iW, $daW, $daVt);
      // Sort the singular values in increasing order
      ReorderSVD($daA, $iH, $iW, $daW, $daVt);

      // Sum each first set of tensor product components
      $iaComponents = array(0, 1, 8, 24, 64);
      for ($iIndex = 0; $iIndex < 4; ++$iIndex) {
        for ($iC = $iaComponents[$iIndex]; $iC < $iaComponents[$iIndex + 1]; ++$iC) {
          for ($j = 0; $j < $iH; ++$j) {
            for ($i = 0; $i < $iW; ++$i) {
              $daAccum[$j*$iW + $i] += $daA[$j*$iW + $iC]*$daW[$iC]*$daVt[$i*$iW + $iC];
            }
          }
          for ($j = 0; $j < $iH; ++$j) {
            for ($i = 0; $i < $iW; ++$i) {
              $iPixel = IntVal(Clamp(Round($daAccum[$j*$iW + $i])));
              ImageSetPixel($qaChannelImages[$iChannel], $i, $j + ($iIndex + 1)*$iH, $iPixel);
            }
          }
        }
      }
      // Make the full color image last.
      for ($j = 0; $j < 5*$iH; ++$j) {
        for ($i = 0; $i < $iW; ++$i) {
          $iRed    = ImageColorAt($qImageR, $i, $j);
          $iGreen  = ImageColorAt($qImageG, $i, $j);
          $iBlue   = ImageColorAt($qImageB, $i, $j);
          $iColor  = (($iRed << 16) + ($iGreen << 8) + $iBlue);
          ImageSetPixel($qImageColor, $i, $j, $iColor);
        }
      }
    }

    $qSvdChannelImage = ImageCreateTrueColor(4*$iW, 5*$iH);
    // Copy the transformed color channel images
    ImageCopyMerge($qSvdChannelImage, $qImageColor, 0, 0, 0, 0, $iW, 5*$iH, 100);
    ImageCopyMerge($qSvdChannelImage, $qImageR, $iW, 0, 0, 0, $iW, 5*$iH, 100);
    ImageCopyMerge($qSvdChannelImage, $qImageG, 2*$iW, 0, 0, 0, $iW, 5*$iH, 100);
    ImageCopyMerge($qSvdChannelImage, $qImageB, 3*$iW, 0, 0, 0, $iW, 5*$iH, 100);

    ImageDestroy($qImageColor);
    ImageDestroy($qImageR);
    ImageDestroy($qImageG);
    ImageDestroy($qImageB);
    return $qSvdChannelImage;
  }
  
  // Use a shell sort to arrange the singular values in increasing order.
  function ReorderSVD($daAtoU, $iM, $iN, $daW, $daVt) {
    $iIncr = 1;
    $daSU = new SplFixedArray($iM);
    $daSV = new SplFixedArray($iN);
    // Find the first increment larger than N
    do {
      $iIncr *= 3;
      ++$iIncr; 
    } while ($iIncr <= $iN);
    do {
      $iIncr = intdiv($iIncr, 3);
      for ($i = $iIncr; $i < $iN; ++$i) {
        $dSW = $daW[$i];
        for ($k = 0; $k < $iM; ++$k) {
          $daSU[$k] = $daAtoU[$k*$iN + $i];
        }
        for ($k = 0; $k < $iN; ++$k) {
          $daSV[$k] = $daVt[$k*$iN + $i];
        }
        $j = $i;
        while ($daW[$j-$iIncr] < $dSW) {
          $daW[$j] = $daW[$j-$iIncr];
          for ($k = 0; $k < $iM; ++$k) {
            $daAtoU[$k*$iN + $j] = $daAtoU[$k*$iN + $j-$iIncr];
          }
          for ($k = 0; $k < $iN; ++$k) {
            $daVt[$k*$iN + $j] = $daVt[$k*$iN + $j-$iIncr];
          }
          $j -= $iIncr;
          if ($j < $iIncr) {
            break;
          }
        }
        $daW[$j] = $dSW;
        for ($k = 0; $k < $iM; ++$k) {
          $daAtoU[$k*$iN + $j] = $daSU[$k];
        }
        for ($k = 0; $k < $iN; ++$k) {
          $daVt[$k*$iN + $j] = $daSV[$k];
        }
      }
    } while ($iIncr > 1);
    for ($k = 0; $k < $iN; ++$k) {
      $s = 0;
      for ($i = 0; $i < $iM; ++$i) {
        if ($daAtoU[$i*$iN + $k] < 0.0) {
          ++$s;
        }
      }
      for ($j = 0; $j < $iN; ++$j) {
        if ($daVt[$j*$iN + $k] < 0.0) {
          ++$s;
        }
      }
      if ($s > ($iM + $iN)/2.0) {
        for ($i = 0; $i < $iM; ++$i) {
          $daAtoU[$i*$iN + $k] = -$daAtoU[$i*$iN + $k];
        }
        for ($j = 0; $j < $iN; ++$j) {
          $daVt[$j*$iN + $k] = -$daVt[$j*$iN + $k];
        }
      }
    }
  }

  $qOriginalImage = ImageCreateFromPng("ChristPantocrator.png");

  $qSvdChannelImage = CreateSvdColorChannelsImage($qOriginalImage);

  Header('Content-type: image/png');
  ImagePng($qSvdChannelImage);
  ImageDestroy($qOriginalImage);
  ImageDestroy($qSvdChannelImage);
?>

Output