function Target_coord_ijk = LDDMM_FindCorrespondingFibersV1p3(Template_LDDMM_dirs, Template_LDDMM_files, Target_LDDMM_dirs, Target_dirs, Target_PreProc_files, working_path,lowerThreshold,upperThreshold)

% function Target_coord_ijk = LDDMM_FindCorrespondingFibersV1p3(Template_LDDMM_dirs, Template_LDDMM_files, Target_LDDMM_dirs, Target_dirs, Target_PreProc_files, working_path,lowerThreshold,upperThreshold)
%
% This function maps a region interest that is labeled on the Atlas
% coordinate to the corresponding location on the target fiber/sheet
% geometry. The region of interest is defined using an intensity treshhold
% range (lowwer to upper). The i,j,k indices of this region of interst is
% then traced back to the i,j,k indices of the target image to be used for
% identifying fiber and sheet angle information.
% 
% 
% 
%  1) Select points in the template anatomy, IJK  at LDDMM match resolution 120x120x110;
%  2) Convert each of these points to points in the Target anatomy, I'J'K' at
%           LDDMM match resolution 120x120x110
%  3) Convert each of these points to points in the Target anatomy, I"J"K" at
%           imaged matched resolution 256x256x117
%
% 
% Inputs:   
% 
%       Template_LDDMM_dirs: cell array representing full path directory that contains template file stored as
%                            xxxx.dat Ex (in linux environment): {'/somePath/Template_It3'};
% 
%       Template_LDDMM_files: cell array listing the template file name Ex.: {'Atlas.dat'};
%
%       Target_LDDMM_dirs: cell array contating directory name for each target. Within these directories. one can find  
%                          LDDMM based transformation files
%                          (Kixmap.dat,...) Ex: 
%                          Target_LDDMM_dirs = {'/somePath/Target_080803';
%                          '/somePath/Target_101003';
%                          '/somePath/Target_101703';
%                          '/somePath/Target_102403';
%                          '/somePath/Target_103103';
%                          '/somePath/Target_112103';
%                          '/somePath/Target_042604';
%                          '/somePath/Target_082303';
%                          '/somePath/Target_112803';
%                          '/somePath/Target_040204';
%                          '/somePath/Target_041804';
%                          '/somePath/Target_050204';};    
% 
% 
%       Target_dirs: cell array containg full path to the target directories
%                 that contain information about rigid body matching and resmpling.
%                 Ex: 
%                 Target_dirs =  {'/somePath/CA_ANALYSIS/Exp_0808';
%                 '/somePath/CA_ANALYSIS/Exp_1010';
%                 '/somepath/CA_ANALYSIS/Exp_1017';
%                 '/somePath/CA_ANALYSIS/Exp_1024';
%                 '/somePath/CA_ANALYSIS/Exp_1031';
%                 '/somePath/CA_ANALYSIS/Exp_1121';
%                 '/somePath/CA_ANALYSIS/Exp_0426';
%                 '/somePath/CA_ANALYSIS/Exp_0823';
%                 '/somePath/CA_ANALYSIS/Exp_1128';
%                 '/somePath/CA_ANALYSIS/Exp_0402';
%                 '/somePath/CA_ANALYSIS/Exp_0418';
%                 '/somePath/CA_ANALYSIS/Exp_0502';};
% 
% 
%       Target_PreProc_files: cell array containing file name that store
%                         preprocessing information such as croppping, resampling, and rigid
%                         transformation.
%                         Ex: 
%                         Target_PreProc_files = {'InvLDDMMSummary_0808.txt';
%                         'InvLDDMMSummary_1010.txt';
%                         'InvLDDMMSummary_1017.txt';
%                         'InvLDDMMSummary_1024.txt';
%                         'InvLDDMMSummary_1031.txt';
%                         'InvLDDMMSummary_1121.txt';
%                         'InvLDDMMSummary_0426.txt';
%                         'InvLDDMMSummary_0823.txt';
%                         'InvLDDMMSummary_1128.txt';
%                         'InvLDDMMSummary_0402.txt';
%                         'InvLDDMMSummary_0418.txt';
%                         'InvLDDMMSummary_0502.txt';};
% 
% 
%        working_path: A string representing the full path for the working
%                directory. This directory will be used to write results on the hard
%                drive, if it is needed.
%                Ex:
%                working_path = '/somePath/LVworkSpace';
%              
%        lower/upperThresholds: Thresholds between 0 to 200 to determine a
%                               region of interest on the template geometry.
% 
%  output: 
%         Target_coord_ijk: cell array, elements of which represent individual target geometry i,j,k indices of a common region of interest 
%                           that is defined on the template geometry.  
%                           
% 
% 
% ICM Oct. 28, 2011




if ~all(~cellfun(@isempty,Template_LDDMM_dirs))
    
    error('Error: Template_LDDMM_dirs information is missing.');
    
end


if ~all(~cellfun(@isempty,Template_LDDMM_files))
    
    error('Error: Template_LDDMM_files information is missing.');
        
end

if ~all(~cellfun(@isempty,Target_LDDMM_dirs))
    
    error('Error: Target_LDDMM_dirs information is missing.');
        
end

if ~all(~cellfun(@isempty,Target_dirs))
    
    error('Error: Target_dirs information is missing.');
        
end

if ~all(~cellfun(@isempty,Target_PreProc_files))
    
    error('Error: Target_PreProc_files information is missing.');
        
end

if isempty(working_path)
    
    error('Error: working path is not provided.');
                 
end


if isempty(lowerThreshold)
    
    error('Error: lower threshold is not provided.');
                 
end


if isempty(upperThreshold)
    
    error('Error: upper threshold is not provided.');
                 
end



if (numel(Target_LDDMM_dirs) ~= numel(Target_dirs)   || numel(Target_LDDMM_dirs) ~= numel(Target_PreProc_files) || numel(Target_PreProc_files) ~= numel(Target_dirs) )
    
    error('Error: number of elements for target LDDMM directories, target pre_processing files, or target directories are not consistent.');
    
end



if ~all(cellfun(@ischar,Template_LDDMM_dirs))
    
    error('Error: Template_LDDMM_dirs information should be provided as string.');
    
end


if ~all(cellfun(@ischar,Template_LDDMM_files))
    
    error('Error: Template_LDDMM_files information should be privided as string.');
        
end

if ~all(cellfun(@ischar,Target_LDDMM_dirs))
    
    error('Error: Target_LDDMM_dirs information should be provided as string.');
        
end

if ~all(cellfun(@ischar,Target_dirs))
    
    error('Error: Target_dirs information should be provided as string.');
        
end

if ~all(cellfun(@ischar,Target_PreProc_files))
    
    error('Error: Target_PreProc_files information should be provided as string.');
        
end



Target_LDDMM_files = {'Kixmap.dat';
    'Kiymap.dat';
    'Kizmap.dat'};



%%%%%%%%%%%%%
% First load the Template Anatomy

fprintf('Loading Template %s %s\n', Template_LDDMM_dirs{1}, Template_LDDMM_files{1});
cd(Template_LDDMM_dirs{1});
Atlas = load_dat_file(Template_LDDMM_files{1},'b');


template_ijk = [];

[template_ijk(:,1), template_ijk(:,2), template_ijk(:,3)] = ind2sub(size(Atlas),find((Atlas > (lowerThreshold)) & (Atlas < (upperThreshold))));

    
    cd(working_path);
    
    
    
    for data_set = 1:length(Target_LDDMM_dirs),
        
        %%%%%%   Step 2 %%%%%
        %  Convert each of these points to points in the Target anatomy, I'J'K' at
        %           LDDMM match resolution 120x120x110
        
        cd(Target_LDDMM_dirs{data_set});
        fprintf('Loading Target %s %s\n', Target_LDDMM_dirs{data_set}, Target_LDDMM_files{1});
        kix = load_dat_file(Target_LDDMM_files{1},'b');
        fprintf('Loading Target %s %s\n', Target_LDDMM_dirs{data_set}, Target_LDDMM_files{2});
        kiy = load_dat_file(Target_LDDMM_files{2},'b');
        fprintf('Loading Target %s %s\n', Target_LDDMM_dirs{data_set}, Target_LDDMM_files{3});
        kiz = load_dat_file(Target_LDDMM_files{3},'b');
        
        template_index = sub2ind(size(Atlas), template_ijk(:,1), template_ijk(:,2), template_ijk(:,3));
        
        fprintf('Interpolating Locations\n'); 
        target_ijk = [kix(template_index)+1, kiy(template_index)+1, kiz(template_index)+1];
        

        

        
        clear kix kiy kiz;
        
        %%%%%%   Step 3 %%%%%
        %   Convert each of these points to points in the Target anatomy, I"J"K" at
        %           imaged matched resolution 256x256x117
        
        fprintf('Loading Target %s %s\n', Target_dirs{data_set}, Target_PreProc_files{data_set});
        cd(Target_dirs{data_set});
        
        fid = fopen(Target_PreProc_files{data_set},'rt');
        done = 0;
        numOperations = 0;
        
        I = target_ijk(:,1); J = target_ijk(:,2); K = target_ijk(:,3);
       
        while(~done)
            op = fgetl(fid);
            if (op == -1)
                done = 1;
            else
                switch(op)
                    case 'Break'
                        done = 1;
                        % fprintf('Break\n');
                    case 'Crop'
                        cropWord = fgetl(fid);
                        cropFactor = str2num(cropWord);    
                        
                        tI = find(I > cropFactor(1));
                        I(tI) = []; J(tI) = []; K(tI) = [];

                        
                        if ~isempty(template_index)
                            template_index(tI) = [];
                        end
                        
                        tI = find(J > cropFactor(2));
                        I(tI) = []; J(tI) = []; K(tI) = [];

                       if ~isempty(template_index)
                            template_index(tI) = [];
                       end
                        tI = find(K > cropFactor(3));
                        I(tI) = []; J(tI) = []; K(tI) = [];
                        
                        if ~isempty(template_index)

                        template_index(tI) = [];
                        end
                        
                        tSetSize = cropFactor;
                         fprintf('   Crop %d %d %d\n', cropFactor); 
                         fprintf('        Temp Size = %d %d %d \n\n', tSetSize);
                        
                    case 'Pad'
                        padWord = fgetl(fid);
                        padFactor = str2num(padWord);    
                        
                        tSetSize = padFactor;
                         fprintf('   Pad %d %d %d\n', padFactor); 
                         fprintf('       Temp Size = %d %d %d \n\n', tSetSize);
                        
                    case 'PadTranslateRoll'
                        padTranslateRollWord = fgetl(fid);
                        padTranslateRollFactor = str2num(padTranslateRollWord);    
                        
                        % Pad translates I,J,K
                        I = I + padTranslateRollFactor(1);
                        J = J + padTranslateRollFactor(2);
                        K = K + padTranslateRollFactor(3);
                        tSetSize = tSetSize + 2*padTranslateRollFactor(1:3);
                        
                        % Translate move coordinates
                        I = I + padTranslateRollFactor(4);
                        J = J + padTranslateRollFactor(5);
                        K = K + padTranslateRollFactor(6);
                        
                        % Roll
                        if (padTranslateRollFactor(7) > 0)
                            tI = find(K > tSetSize(3));
                            K(tI) = K(tI)- tSetSize(3); 
                        else
                            tI = find( K < 0.5);
                            K(tI) = tSetSize(3) - abs(K(tI)); 
                        end
                        
                          fprintf('   Pad/Translate/Roll %d %d %d %d %d %d %d\n', padTranslateRollFactor);
                          fprintf('        Temp Size = %d %d %d \n\n', tSetSize);
                    case 'Roll'
                        rollWord = fgetl(fid);
                        rollFactor = str2num(rollWord); 
                        
                        K = K + rollFactor;
                        
                        if (rollFactor > 0)
                            tI = find(K > tSetSize(3));
                            K(tI) = K(tI)- tSetSize(3); 
                        else
                            tI = find( K < 0.5);
                            K(tI) = tSetSize(3) - abs(K(tI)); 
                        end
                    case 'Resize'
                        resizeWord = fgetl(fid);
                        resizeFactor = str2num(resizeWord);    
                        
                        I = I * resizeFactor(1);
                        J = J * resizeFactor(2);
                        K = K * resizeFactor(3);
                        
                        tSetSize = tSetSize .* resizeFactor;
                         fprintf('   Resize %0.4f %0.4f %0.4f\n', resizeFactor); 
                         fprintf('        Temp Size = %d %d %d \n\n', tSetSize);
                    case 'Translate'
                        translateWord = fgetl(fid);
                        translateFactor = str2num(translateWord);   
                        
                        % Translate move coordinates
                        I = I + translateFactor(1);
                        J = J + translateFactor(2);
                        K = K + translateFactor(3);
                        
                        fprintf('   Translate %d %d %d\n\n', translateFactor); 
                        
                    case 'RigidTransform'
                        rotationFileName = fgetl(fid);
                        translationFileName = fgetl(fid);
                        R = load(rotationFileName);
                        T = load(translationFileName);
                        [IJK] = (R' * ([I(:)-1,J(:)-1,K(:)-1]' - repmat(T(:),1,length(I))))';
                        
                        I = IJK(:,1)+1;
                        J = IJK(:,2)+1;
                        K = IJK(:,3)+1;
                        clear IJK;
                        
                         fprintf('   Rigid Transform %s %s\n\n', rotationFileName, translationFileName);
                        
                    case 'InvRigidTransform'
                        rotationFileName = fgetl(fid);
                        translationFileName = fgetl(fid);
                        R = load(rotationFileName);
                        T = load(translationFileName);
                        [IJK] = (R * [I(:)-1,J(:)-1,K(:)-1]' + repmat(T(:),1,length(I)))';
                        
                        I = IJK(:,1)+1;
                        J = IJK(:,2)+1;
                        K = IJK(:,3)+1;
                        clear IJK;
                        
                         fprintf('   Inv Rigid Transform %s %s\n\n', rotationFileName, translationFileName);
                        
                    case 'Flip'
                        t = I;
                        I = J;
                        J = t;
                        K = tSetSize(3)-K;
                end
                numOperations = numOperations+1;
                
                if (numOperations > 60) 
                    done = 1; 
                end;
            end
        end
        fclose(fid);
        target_ijk = [I,J,K];
        Target_coord_ijk{data_set} = target_ijk;
        
        

        
        target_ijk = [];
        
        fprintf('-------------------- Heart  Complete -----------------------\n');
       
    
    end
    
    


cd(working_path)




