nonsubsampled contourlet transform for color images

I am fusing panchromatic and multispectral images. I am using nonsubsampled contourlet transform for decomposing the image but NSCT is decomposing only grayscale image what modification is to be done in the code to decompose the color image. please heip me out in this.

iam using the following code.

clc;
close all;
clear all;
path(path,'PCNN_toolbox/')
path(path,'nsct_toolbox')
path(path,'FusionEvaluation/')
Low_Coeffs_Rule='SF-PCNN'; %'ave','PCNN','SF-PCNN'
High_Coeffs_Rule='SF-PCNN'; %'max','PCNN','SF-PCNN'
ori_A=imread('pan.jpg');
ori_A = rgb2gray(ori_A);
ori_B=imread('ms.jpg');
ori_B = rgb2gray(ori_B);
A=double(ori_A)/255;
B=double(ori_B)/255;
%%Parameters for NSCT
pfilt = '9-7';
dfilt = 'pkva';
nlevs = [0,1,3,4,4];% default
%%Parameters for PCNN
Para.iterTimes=200;
Para.link_arrange=3;
Para.alpha_L=1;% 0.06931 Or 1
Para.alpha_Theta=0.2;
Para.beta=3;% 0.2 or 3
Para.vL=1.0;
Para.vTheta=20;
disp('Decompose the image via nsct ...')
yA=nsctdec(A,nlevs,dfilt,pfilt);
yB=nsctdec(B,nlevs,dfilt,pfilt);
%%%%%%%%%%%%save the NSCT coefficients %%%%%%%%%%%%
 save yA yA
 save yB yB
%%%%%%%%%%%%load the NSCT coefficients %%%%%%%%%%%%
load yA
load yB
n = length(yA);
%%Initialized the coefficients of fused image
Fused=yA;
%=============================================
% Lowpass subband
disp('Process in Lowpass subband...')
ALow1= yA{1};
BLow1 =yB{1};
ALow2= yA{2};
BLow2 =yB{2};
switch Low_Coeffs_Rule
    case 'average'
        Fused{1}=(ALow1+BLow1)/2;
        Fused{2}=(ALow2+BLow2)/2;
    case 'PCNN'       
        Fused{1}=fusion_NSCT_PCNN(ALow1,BLow1,Para);
        Fused{2}=fusion_NSCT_PCNN(ALow2,BLow2,Para);
    case 'SF-PCNN'
       Fused{1}=fusion_NSCT_SF_PCNN(ALow1,BLow1,Para);
       Fused{2}=fusion_NSCT_SF_PCNN(ALow2,BLow2,Para);
end
%=============================================
% Bandpass subbands
disp('Process in  Bandpass subbands...')
for l = 3:n
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for d = 1:length(yA{l})
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Ahigh = yA{l}{d};
        Bhigh = yB{l}{d};
        switch High_Coeffs_Rule
            case 'max'
                decision_map=(abs(Ahigh)>=abs(Bhigh));
                Fused{l}{d}=decision_map.*Ahigh + (~decision_map).*Bhigh;
            case 'PCNN'
                Fused{l}{d}=fusion_NSCT_PCNN(Ahigh,Bhigh,Para);
            case 'SF-PCNN'
                Fused{l}{d}=fusion_NSCT_SF_PCNN(Ahigh,Bhigh,Para);
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end
end
disp('High frequecy field process is ended')
%=============================================
disp('Reconstruct the image via nsct ...')
F=nsctrec(Fused, dfilt, pfilt);
disp('Reconstruct is ended...')
F=F*255;
F(F<0)=0;
disp('F>255')
F(F>255)=255;
F=round(F);
Eval=Evaluation(double(ori_A),double(ori_B),F,256)
figure,imshow(ori_A,[])
figure,imshow(ori_B,[])
figure,imshow(F,[])
figure,imshow(F-double(ori_A),[])
figure,imshow(F-double(ori_B),[])

[angel]

hey... m working on the same field... did u find d solution? i am trying in different set of algorithms but i am not getting any solutions.