function [bestWaveletsAll, PCANet] = findBestWaveletsTesting(imagesCell, orient_default, orient_best, PCANet, numImagesToUse, imageSize, fidLogs, param, stepPrint, plotta)

%init

bestWaveletsAll(param.numBestWavelets).filter = [];

%display

fprintf_pers(fidLogs, ['\t\t\t' num2str(numImagesToUse) ' images used for Gabor analysis... \n'])

%--------------------------------------

%Compute parametrized wavelets

%parametrizedWavelets = gaborArrayParametrized(param.sigma, param.wavelength, param.divThetaParametrized);

parametrizedWavelets = gaborArrayParametrized(param.sigma, param.wavelength, deg2rad(orient_default));

%--------------------------------------

%Compute complete multi-scale Gabor filter bank

%fprintf_pers(fidLogs, '\t\tComputing Gabor array... \n')

gaborBank = gaborArrayFromScales(imageSize, deg2rad(unique([orient_default orient_best])), param, plotta);

%init counter for how many times each wavelet is chosen

o_counterAll = zeros(numel(gaborBank), 1);

    if plotta

        fsfigure

        for indf = 1 : param.divThetaParametrized

            %sum(real(parametrizedWavelets(indf).filter(:)))

            subplot(1,param.divThetaParametrized,indf)

            imshow(real(parametrizedWavelets(indf).filter),[])

            axis on

            axis image

        end %for indf

        suplabel([num2str(param.divThetaParametrized) ' filters fixed']);

        figure,

        imhist(real(parametrizedWavelets(indf).filter))

        %pause

    end %if plotta

%loop on images

parfor j = 1 : numImagesToUse

% for j = 1 : numel(imagesCell)

    %get id of current worker

    t = getCurrentTask();

    %display progress

    if mod(j, stepPrint) == 0

        fprintf(1, ['\t\tCore ' num2str(t.ID) ': ' num2str(j) ' / ' num2str(numImagesToUse) '\n'])

    end %if mod(i, 100) == 0

    %get img

    im = imagesCell{j};

    %test img alternative

    %%%%%%%%%%%%%%%%%%%%%%%%%%

    %im = loadPreprocessImage('D:\UoT\Projects\Palmprint Gabor tuning\Code\alt\GaborWavelet\images\1_1.jpg', imageSize);

    %im = loadPreprocessImage('D:\UoT\Projects\Palmprint Gabor tuning\Code\alt\GaborWavelet\images\test image\Peppers.tiff', imageSize);

    %%%%%%%%%%%%%%%%%%%%%%%%%%

    %pause

    %--------------------------------------

    %Compute frequency of Gabor filter bank

    %fprintf_pers(fidLogs, '\t\tComputing spectrum of gabor array... \n')

    %totFreq = computeFrequencyGaborArray(gaborBank, plotta);

    %--------------------------------------

    %Perform reference Gabor filtering

    %fprintf_pers(fidLogs, '\t\tPerforming reference Gabor filtering... \n')

    resF = referenceGaborFilter(gaborBank, im);

    %--------------------------------------

    %create powerMap for all filters

    %fprintf_pers(fidLogs, '\t\tComputing powermaps... \n')

    powerMapF = computePowerMaps(resF, gaborBank);

    %--------------------------------------

    %Sort filter response information

    %fprintf_pers(fidLogs, '\t\tSorting filter response information... \n')

    [sortRes, sizeRes] = sortFilterResponse(gaborBank, powerMapF);

    %--------------------------------------

    %Compute most used wavelets

    %fprintf_pers(fidLogs, '\t\t\tComputing most used wavelets... \n')

    o_counter = getMostUsedWavelets(sortRes, gaborBank, powerMapF, param, sizeRes, plotta);

    %increment

    o_counterAll = o_counterAll + o_counter;

end %for j = 1 : length(vectorIndexTest)

%display most used filters

if plotta

    figure,

    bar(1:size(o_counterAll), o_counterAll)

    xlabel('Wavelet n.')

    ylabel('Perc. of occurence')

    title('Most used wavelets');

end %if plotta

%--------------------------------------

%sort most used filters counter

[~, ind_o_counter_All_sort] = sort(o_counterAll, 'descend');

%consider only best wavelets

ind_o_counter_All_sort = ind_o_counter_All_sort(1 : param.numBestWavelets);

bestWavelets = gaborBank(ind_o_counter_All_sort);

%display bestWavelets

if plotta

    fsfigure

    for indf = 1 : numel(bestWavelets)

        subplot(1, numel(bestWavelets),indf)

        imshow(bestWavelets(indf).even,[])

        axis on

        axis image

    end %for indf

    suplabel([num2str(numel(bestWavelets)) ' filters most used']);

end %if plotta

%--------------------------------------

%assign to global structure

countfilter = 1;

%parametrized

for o = 1 : numel(parametrizedWavelets)

    bestWaveletsAll(countfilter).filter = parametrizedWavelets(o).filter;

    countfilter = countfilter + 1;

end %for o

%dynamic

for o = 1 : numel(bestWavelets)

    bestWaveletsAll(countfilter).filter = complex(bestWavelets(o).even, bestWavelets(o).odd);

    countfilter = countfilter + 1;

end %for o

%display

%     if plotta

%         bestWaveletsCell = squeeze(struct2cell(bestWavelets));

%         figure,

%         subplot(1,2,1)

%         montage(bestWaveletsCell(1,:), 'DisplayRange', [])

%         axis on

%         axis image

%         title('Best wavelets: even');

%         subplot(1,2,2)

%         montage(bestWaveletsCell(2,:), 'DisplayRange', [])

%         title('Best wavelets: odd');

%         pause

%     end %if plotta

%--------------------------------------

%loop on images, test reconstruct

if 0

    for j = 1 : numel(imagesCell)

        %get img

        im = imagesCell{j};

        %Perform reference Gabor filtering

        resF = referenceGaborFilter(bestWaveletsAll, im);

        %create powerMap for all filters

        powerMapF = computePowerMaps(resF, bestWaveletsAll);

        %Sort filter response information

        [sortRes, sizeRes] = sortFilterResponse(bestWaveletsAll, powerMapF);

        %--------------------------------------

        %Loop on wavelets - all filters

        %fprintf_pers(fidLogs, '\t\t\tReconstructing image with all filters... \n')

        % allFiltersIndexes = 1:numel(gaborBank);

        %[tmpImage, errV] = reconstructImage(resF, sortRes, gaborBank, powerMapF, allFiltersIndexes, im, imageSize, sizeRes, param, plotta);

        %--------------------------------------

        %Loop on wavelets - best filters

        %fprintf_pers(fidLogs, ['\t\t\tReconstructing image with the best ' num2str(param.numBestWavelets) ' filters... \n'])

        [tmpImageBest, errBestV] = reconstructImage(resF, sortRes, bestWaveletsAll, powerMapF, [1 : numel(bestWaveletsAll)], im, imageSize, sizeRes, param, plotta);

        pause

    end %for j = 1 : numel(imagesCell)

end %if 0

%--------------------------------------

%update number of filters

%assign updated number of filters to PCANet

PCANet.NumFilters(end) = numel(bestWaveletsAll);