Decompose polarimetric SAR covariance matrices in MATLAB

I have the following code to implement the algorithm described in the article [_Adaptive Model-Based Decomposition of Polarimetric SAR Covariance Matrices_](<http://dx.doi.org/10.1109/TGRS.2010.2076285>) by Arii, van Zyl, and Kim, in IEEE Transactions on Geoscience and Remote Sensing.

I want to apply it to 3 images of the sizes 676×718, 1430×1440, and 5576×9444, but the code is too slow even for the 676×718 image.

    function [Ps,Pd,Pv] = Arii2011_Modified(C11,C12_imag,C12_real,C13_imag,C13_real,C22,C23_imag,C23_real,C33)
    [Ps,Pd,Pv] = arrayfun(@Arii2011_Modified_1Pixel,C11,C12_imag,C12_real,C13_imag,C13_real,C22,C23_imag,C23_real,C33);
    end

--

    function [Ps,Pd,Pv] = Arii2011_Modified_1Pixel(C11,C12_imag,C12_real,C13_imag,C13_real,C22,C23_imag,C23_real,C33)
    global DataInstance;
    DataInstance = Data(C11,C12_imag,C12_real,C13_imag,C13_real,C22,C23_imag,C23_real,C33);
    MeanOrientationAngleStep = 1/100;
    StandardDeviationStep = 1/100;
    NumberOfAngles = floor(2*pi/MeanOrientationAngleStep);
    NumberOfStandandardDeviations = 1/StandardDeviationStep;
    theta = -pi+MeanOrientationAngleStep*(0:NumberOfAngles)';
    sigma = StandardDeviationStep*(0:NumberOfStandandardDeviations)';
    Pv = FindFv(sigma,theta);
    Ps = 0;
    Pd = 0;
    end

--

    function Fv = FindFv(sigmav,thetav)
    global DataInstance;
    [sigma,theta] = meshgrid(sigmav,thetav);
    theta = reshape(theta,[],1);
    sigma = reshape(sigma,[],1);
    coef1 = 1/2-((((1-2*sigma.^2).*(1-sigma.^2).*cos(2*theta))./(2*(1+sigma.^2))));
    coef2 = ((1-4*sigma.^2).*(1-3*sigma.^2).*(1-2*sigma.^2).*(1-sigma.^2).*cos(4*theta))./(8*(1+sigma.^2).*(1+2*sigma.^2).*(1+3*sigma.^2));
    coef3 = ((1-2*sigma.^2).*(1-sigma.^2).*sin(2*theta))./(2*sqrt(2)*(1+sigma.^2));
    coef4 = ((1-4*sigma.^2).*(1-3*sigma.^2).*(1-2*sigma.^2).*(1-sigma.^2).*sin(4*theta))./(4*sqrt(2)*(1+sigma.^2).*(1+2*sigma.^2).*(1+3*sigma.^2));
    Fv_calc(:,1) = DataInstance.C11./(coef1+coef2);
    aQuadratic1 = -2*coef2.*(coef1+coef2);
    aQuadratic2 = (coef3-coef4).^2;
    aQuadratic = aQuadratic1 - aQuadratic2;
    bQuadratic1 = 2*DataInstance.C11*coef2-DataInstance.C22*(coef1+coef2);
    bQuadratic2 = -2*DataInstance.C12_real*(coef3-coef4);
    bQuadratic = bQuadratic1 - bQuadratic2;
    cQuadratic1 = DataInstance.C11*DataInstance.C22;
    cQuadratic2 = DataInstance.C12_real^2+DataInstance.C12_imag^2;
    cQuadratic = cQuadratic1 - cQuadratic2;
    for i = 1 : size(sigma,1)
        rQuadratic = roots([aQuadratic(i,1) bQuadratic(i,1) cQuadratic]);
        sel1 = rQuadratic == real(rQuadratic);
        rQuadratic = rQuadratic(sel1);
        sel2 = rQuadratic == abs(rQuadratic);
        rQuadratic = rQuadratic(sel2);
        Fv_calc(i,2) = max(rQuadratic);
    end
    aCubic1 = 2*coef2.*(1-coef1+coef2).*(coef1+coef2);
    aCubic2 = 2*coef2.*(coef3-coef4).*(coef3+coef4);
    aCubic3 = -(coef1+coef2).*(coef3+coef4).^2;
    aCubic4 = 2*coef2.^3;
    aCubic5 = -(1-coef1+coef2).*(coef3-coef4).^2;
    aCubic = aCubic1+aCubic2-aCubic3-aCubic4-aCubic5;
    bCubic1 = -2*coef2.*(DataInstance.C11*(1-coef1+coef2)+DataInstance.C33*(coef1+coef2))+DataInstance.C22*(1-coef1+coef2).*(coef1+coef2);
    bCubic2 = -2*coef2.*(DataInstance.C23_real*(coef3-coef4)+DataInstance.C12_real*(coef3+coef4))+2*DataInstance.C13_real*(coef3-coef4).*(coef3+coef4);
    bCubic3 = 2*DataInstance.C23_real*(coef1+coef2).*(coef3+coef4)+DataInstance.C11*(coef3+coef4).^2;
    bCubic4 = (4*DataInstance.C13_real+DataInstance.C22)*coef2.^2;
    bCubic5 = 2*DataInstance.C12_real*(1-coef1+coef2).*(coef3-coef4)+DataInstance.C33*(coef3-coef4).^2;
    bCubic = bCubic1+bCubic2-bCubic3-bCubic4-bCubic5;
    cCubic1 = 2*DataInstance.C11*DataInstance.C33*coef2-DataInstance.C11*DataInstance.C22*(1-coef1+coef2)-DataInstance.C22*DataInstance.C33*(coef1+coef2);
    cCubic2 = 2*(DataInstance.C12_real*DataInstance.C23_real-DataInstance.C12_imag*DataInstance.C23_imag)*coef2-...
        2*(DataInstance.C13_real*DataInstance.C23_real+DataInstance.C13_imag*DataInstance.C23_imag)*(coef3-coef4)-...
        2*(DataInstance.C12_real*DataInstance.C13_real+DataInstance.C12_imag*DataInstance.C13_imag)*(coef3+coef4);
    cCubic3 = -(DataInstance.C23_real^2+DataInstance.C23_imag^2)*(coef1+coef2)-2*DataInstance.C11*DataInstance.C23_real*(coef3+coef4);
    cCubic4 = 2*(DataInstance.C13_real^2+DataInstance.C13_imag^2+DataInstance.C13_real*DataInstance.C22)*coef2;
    cCubic5 = -(DataInstance.C12_real^2+DataInstance.C12_imag^2)*(1-coef1+coef2)-2*DataInstance.C33*DataInstance.C12_real*(coef3-coef4);
    cCubic = cCubic1+cCubic2-cCubic3-cCubic4-cCubic5;
    dCubic1 = DataInstance.C11*DataInstance.C22*DataInstance.C33;
    dCubic2 = 2*(DataInstance.C12_real*(DataInstance.C13_real*DataInstance.C23_real+DataInstance.C13_imag*DataInstance.C23_imag)+...
        DataInstance.C12_imag*(DataInstance.C13_imag*DataInstance.C23_real-DataInstance.C13_real*DataInstance.C23_imag));
    dCubic3 = DataInstance.C11*(DataInstance.C23_real^2+DataInstance.C23_imag^2);
    dCubic4 = DataInstance.C22*(DataInstance.C13_real^2+DataInstance.C13_imag^2);
    dCubic5 = DataInstance.C33*(DataInstance.C12_real^2+DataInstance.C12_imag^2);
    dCubic = dCubic1+dCubic2-dCubic3-dCubic4-dCubic5;
    for i = 1 : size(sigma,1)
        rCubic = roots([aCubic(i,1) bCubic(i,1) cCubic(i,1) dCubic]);
        sel3 = rCubic == real(rCubic);
        rCubic = rCubic(sel3);
        sel4 = rCubic == abs(rCubic);
        rCubic = rCubic(sel4);
        Fv_calc(i,3) = max(rCubic);
    end
    Fv = min(Fv_calc,[],2);
    Fv = max(Fv);
    end  
Where `DataInstance` is a class as follows:  

    classdef Data < handle
        properties (SetAccess = public)
        C11;
        C12_imag;
        C12_real;
        C13_imag;
        C13_real;
        C22;
        C23_imag;
        C23_real;
        C33;    
        end
        methods (Access = public)
            function obj = Data(C11,C12_imag,C12_real,C13_imag,C13_real,C22,C23_imag,C23_real,C33)
                obj.C11 = C11;
                obj.C12_imag = C12_imag;
                obj.C12_real = C12_real;
                obj.C13_imag = C13_imag;
                obj.C13_real = C13_real;
                obj.C22 = C22;
                obj.C23_imag = C23_imag;
                obj.C23_real = C23_real;
                obj.C33 = C33;
            end
        end
    end  
Here is the result of profiling the code for one pixel:  
[![Profile screenshot][1]][1]

And here is the result when I've tried to run the code for the whole 676×718 image but because of taking a lot of time, I had to interrupt the running process. [![Profiler screenshot][2]][2]

From the first screenshot, it seems that the most problematic part of `FindFv` function is finding the roots but the second image shows that the other arithmetic operations in `FindFv` are a major problem, too.

 - Which is better? Write the whole `FindFv` function as a MEX file or just the root parts?

 - Will writing as a MEX file help a lot considering that I can't use most of those vectorized codes in other parts of `FindFv` other than `roots`?

 - Is there anyway to accelerate the process through MEX files or another way? Will this help significantly?

[1]: //i.stack.imgur.com/IZkbQ.png
[2]: //i.stack.imgur.com/qoxTu.png
[3]: http://i.stack.imgur.com/8807G.png
[4]: http://i.stack.imgur.com/DNpJb.png

[ROY THANKACHAN]

How do you obtained the Polarimetric SAR Covariance Matrices ?

[Juanping Zhao]

I am confusing how to obtain the Polarimetric SAR Coherence Matrice? Nearlly the same question for you. Can you tell me how to solve this problem? Very thanks for you.

Answer 1

[najoua elhajjaji]

i have the same question, could you plz tell me how to find the PolSar covariance matrix using matlab?? Thanks in advance for your help.