Image Processing Computer Science

Singular Value Decomposition on Color Channels

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

      // Sum each first set of tensor product components
      $iaComponents = array(0, 1, 4, 16, 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;
  }

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

  $qSvdChannelImage = CreateSvdColorChannelsImage($qOriginalImage);

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

Output

 
 

© 2007–2022 XoaX.net LLC. All rights reserved.