5
$\begingroup$

I am working on segmenting a textured image using Gabor Filters. This is based on this paper. There exist a Matlab implementation of this as well. Earlier I managed to do up to this. With the help of this now I come up with the following:

(* Design an array of Gabor Filters *)

{numRows, numCols} = ImageDimensions[img];

wavelengthMin = 4/Sqrt[2];
wavelengthMax = Sqrt[nrows^2 + ncols^2];
n = Floor[Log2[wavelengthMax/wavelengthMin]];
wavelength = Table[{2, 2} 2^i, {i, 0, n - 2}];
deltaTheta = Pi/4;
orientation = Table[i, {i, 0, Pi - deltaTheta, deltaTheta}];

gabormag = 
  Flatten[Table[
    ImageAdjust[
     GaborFilter[img, 1, wavelength[[i]], orientation[[j]]]], {i, 1, 
     Length@wavelength}, {j, 1, Length@orientation}]];
gaborWavelength = 
  Flatten[ConstantArray[
    Table[Norm[wavelength[[i]]], {i, 1, Length@wavelength}], 
    Length@orientation]];

(* Post-processing of Gabor Magnitude images into Gabor Features *)

K = 3;
gabormagfiltered = gabormag;
Table[
  sigma = 0.5 gaborWavelength[[i]];
  gabormagfiltered[[i]] = 
   ImageData[GaussianFilter[gabormag[[i]], K sigma]]
  , {i, 1, Length@gabormag}];

(* Add a map of spatial location information *)

meshgrid[x_List, y_List] := {ConstantArray[x, Length[y]], 
  Transpose@ConstantArray[y, Length[x]]}

{X, Y} = meshgrid[Range[ncols], Range[nrows]];
featureset = Append[Append[gabormagfiltered, X], Y];

(* Reshape Data *)

npoints = nrows*ncols;
freq = Dimensions[featureset][[1]];

Needs["Developer`"]
X = Developer`ToPackedArray[
   ArrayReshape[featureset, {freq, npoints}]];

(* Normalize the features to be zero mean, unit variance *)

Xstandardized = Standardize[X];

(* The following is based on the answer provided by @Michael E2 here) ( Apply Principal Component Analysis *)

pca = PrincipalComponents[Xstandardized];
X1 = Transpose[Transpose@Xstandardized - Mean[Xstandardized]];

(* Extract coefficients *)

{u, \[Sigma], v} = SingularValueDecomposition[X1];

However, when I apply SingularValueDecomposition on X1 my machine (intel core i7, 16 GB RAM) runs out of memory, while in Matlab it works absolutely fine.

How can I efficiently apply PCA and SVD to extract the coefficients that I get directly using pca in Matlab?

For example, below are the original image, feature2DImage and clustered image (with 5 clusters)

enter image description here

enter image description here

enter image description here

$\endgroup$
  • $\begingroup$ If you don't need all of them, try {u, \[Sigma], v} = SingularValueDecomposition[X1, Length@SingularValueList[X1]]; $\endgroup$ – Michael E2 May 8 '18 at 12:41
  • $\begingroup$ @MichaelE2 After getting v following your approach, I applied feature2DImage = ImageAdjust@ImageResize[Image[X1.v], {225, 225}]. However, this does not give the desired image. Where am I doing wrong? $\endgroup$ – Majis May 8 '18 at 12:49
  • $\begingroup$ @MichaelE2 Probably Dimensions[v] = {50625, 21} is a problem. $\endgroup$ – Majis May 8 '18 at 12:51
  • $\begingroup$ I think what you want is PadRight[u.\[Sigma], {Automatic, Length@v}], since for the full SVD, u.\[Sigma] == X1.v. But it still doesn't look like the right image. If you have MATLAB, you might check the intermediate results. Somehow I doubt X1 is supposed to have a dimension 50625 if it's used to reconstruct the image as shown. -- Note Standardize zeros the mean, so Transpose[Transpose@Xstandardized - Mean[Xstandardized]] is unnecessary. $\endgroup$ – Michael E2 May 8 '18 at 13:34
  • 1
    $\begingroup$ Try SingularValueDecomposition[X1, Min@Dimensions@X1] and u.\[Sigma] for coeff, then. -- It would help if there were an example image that everybody could do the same computations as each other. I don't think the web image of the puppy is exactly the same, or maybe you're using a different image. $\endgroup$ – Michael E2 May 8 '18 at 14:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.