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–2024 XoaX.net LLC. All rights reserved.