function [neurites_combined,all_points,run_length] = TrakEM2_XML_extractor(filename,do_smoothing,do_output)

% Extract  relevant information on polylines from XML files generated by
% TrakEM2.  Convert to readable format (3xn matrix of points in the line).

input_filename = filename;

gunzip([input_filename '.xml.gz']);
input_xml = xml_read([input_filename '.xml']);

 img_height = str2num(input_xml.t2_layer_set.t2_layer(1).t2_patch.ATTRIBUTE.height);
 img_width = str2num(input_xml.t2_layer_set.t2_layer(1).t2_patch.ATTRIBUTE.width);
 num_layers = size(input_xml.t2_layer_set.t2_layer,1);
 num_lines = size(input_xml.t2_layer_set.t2_polyline,1)
 
 counters = zeros(num_lines,1);
 
 neurite_stacks = zeros(num_lines,img_height,img_width,num_layers);
 
 
 for i = 1:num_layers
 layer_id{i} = input_xml.t2_layer_set.t2_layer(i).ATTRIBUTE.oid;
 layer_z{i} = input_xml.t2_layer_set.t2_layer(i).ATTRIBUTE.z;
 end
 
run_length = 0;
 
 for k = 1:num_lines
     
     if isfield(input_xml.t2_layer_set.t2_polyline(k).ATTRIBUTE,'d')
 
 line_points = input_xml.t2_layer_set.t2_polyline(k).ATTRIBUTE.d;
 layer_numbers = input_xml.t2_layer_set.t2_polyline(k).ATTRIBUTE.layer_ids;
 transform_matrix = input_xml.t2_layer_set.t2_polyline(k).ATTRIBUTE.transform;
 
 parsed_str_xy = strsplit(line_points);
 parsed_str_z = strsplit(layer_numbers,',');
 str_size = size(parsed_str_xy,2);
 
 parsed_transform_matrix = strsplit(transform_matrix,{',','(',')'});
 x_shift = str2num(parsed_transform_matrix{6});
 y_shift = str2num(parsed_transform_matrix{7});
 
 % Would possibly be smart to get out the z thickness, although you maybe
 % can assume you're looking at 500nm slices (the norm).
 
 coords = zeros(3,str_size/2);
 
for i = 1:(str_size/2)
    parsed_str_xy_2 = strsplit(parsed_str_xy{2*i},',');
    % xcoordinate
    coords(1,i) = str2num(parsed_str_xy_2{1}) + x_shift;
    % y coordinate
    coords(2,i) = str2num(parsed_str_xy_2{2}) + y_shift;
    % z coordinate
    index_mat = strcmp(parsed_str_z{i},layer_id);
    coords(3,i) = find(index_mat == 1);
   % zpoint(i) = str2num(layer_z{index})
end


for i = 1:((str_size/2) - 1)
    run_length = run_length + norm(coords(:,i+1) - coords(:,i));
end
    
    
% generate lines connecting points

% find a run (in x, y, or z) which isn't 0,
% then trace the rise of the line against that run.

for i = 1: ((str_size/2) - 1)
    run = coords(1,i+1) - coords(1,i);
    run2 = coords(2,i+1) - coords(2,i);
    run3 = coords(3,i+1) - coords(3,i);
    
    if run ~= 0
    slope(1) = (coords(2,i+1) - coords(2,i)) / run;
    slope(2) = (coords(3,i+1) - coords(3,i)) / run;
    
    % add in endpoint
    neurite_stacks(k,coords(2,i),coords(1,i),coords(3,i)) = 1;
    counters(k) = counters(k) +1;
    all_points{k}(counters(k),:) = [coords(1,i) coords(2,i) coords(3,i)];
    
    for j = 1:abs(run)
        
        shift = j*sign(run);
        new_x = coords(1,i) + shift;
        new_x_rounded = round(new_x);
        new_y = coords(2,i) + shift*slope(1);
        new_y_rounded = round(new_y);
        new_z = coords(3,i) + shift*slope(2);
        new_z_rounded = round(new_z);
        
        neurite_stacks(k,new_y_rounded,new_x_rounded,new_z_rounded) = 1;
        counters(k) = counters(k) +1;
        all_points{k}(counters(k),:) = [new_x_rounded new_y_rounded new_z_rounded];
    end
    
    end
    
    if run2 ~= 0
        
        slope(1) = (coords(1,i+1) - coords(1,i)) / run2;
        slope(2) = (coords(3,i+1) - coords(3,i)) / run2;
        
        for j = 1:abs(run2)
            
            shift = j*sign(run2);
            new_x = coords(1,i) + shift*slope(1);
            new_x_rounded = round(new_x);
            new_y = coords(2,i) + shift;
            new_y_rounded = round(new_y);
            new_z = coords(3,i) + shift*slope(2);
            new_z_rounded = round(new_z);
            
            neurite_stacks(k,new_y_rounded,new_x_rounded,new_z_rounded) = 1;
            counters(k) = counters(k) +1;
            all_points{k}(counters(k),:) = [new_x_rounded new_y_rounded new_z_rounded];
        end
    end
        
    if run3 ~= 0
        slope(1) = (coords(1,i+1) - coords(1,i)) / run3;
        slope(2) = (coords(2,i+1) - coords(2,i)) / run3;
        
        for j = 1:abs(run3)
            
            shift = j*sign(run3);
            new_x = coords(1,i) + shift*slope(1);
            new_x_rounded = round(new_x);
            new_y = coords(2,i) + shift*slope(2);
            new_y_rounded = round(new_y);
            new_z = coords(3,i) + shift;
            new_z_rounded = round(new_z);
            
            neurite_stacks(k,new_y_rounded,new_x_rounded,new_z_rounded) = 1;
            counters(k) = counters(k) +1;
            all_points{k}(counters(k),:) = [new_x_rounded new_y_rounded new_z_rounded];
        end
    end
            
end
    
     end
end


% convert matrix to maximum for 8 bit images (for use in 3D viewer).
% eight_bit_neurite = uint8(255*neurite_stacks);

brightened_neurites = uint16(70000*neurite_stacks);

%spread out neurites so that you can actually see them to some degree.

if do_smoothing
    neurites_combined = smooth3(squeeze(max(brightened_neurites,[],1)),'gaussian',5,4);
%     stack_max = max(max(max(neurites_combined)));
%     neurites_combined = neurites_combined*255/stack_max;

    neurites_combined = uint16(neurites_combined);
    % neurites_brighter = uint8((neurites>5)*255);
else
    neurites_combined = squeeze(max(brightened_neurites,[],1));
end


if do_output
    TiffStackWrite(neurites_combined,[input_filename],16)
end

end