diff --git a/BIOP_Clustering.m b/BIOP_Clustering.m
new file mode 100644
index 0000000..f509dff
--- /dev/null
+++ b/BIOP_Clustering.m
@@ -0,0 +1,760 @@
+function BIOP_Clustering
+%Start EasyXT
+addpath('EasyXT');
+addpath('functions');
+global X;
+X = EasyXT();
+
+f = figure('Name', 'BIOP - Clustering WingJ', ...
+ 'NumberTitle', 'off', 'Visible', 'off', 'Position', [0 0 700 512], ...
+ 'Tag', 'clustetring_main_f', ...
+ 'DockControls', 'off', 'ToolBar', 'none', 'MenuBar', 'none');
+
+% Move the window to the center of the screen.
+movegui(f,'center');
+
+% Make the window visible.
+f.Visible = 'on';
+
+% Add 1 panel on the left and one on the right
+settingsp = uipanel(f,'Title','Settings', 'Tag', 'p_settings',...
+ 'Position',[.05 .05 .65 .9]);
+
+commandsp = uipanel(f,'Title','Actions', 'Tag', 'p_commands',...
+ 'Position',[.70 .05 .25 .9]);
+
+% Settings Tab
+uitabgroup('Parent', settingsp, 'Tag', 'tab_channels');
+
+% Add buttons for refreshing data, and batch
+h = .1;
+w = .9;
+l = .05;
+
+uicontrol(commandsp,'Style','pushbutton','String','Refresh',...
+ 'Units','normalized', 'Tag', 'btn_refresh', ...
+ 'Position',[l .90 w h],...
+ 'Callback', @btn_refresh_callback);
+uicontrol(commandsp,'Style','pushbutton','String','Get Sum Channel',...
+ 'Units','normalized', 'Tag', 'btn_sum', ...
+ 'Position',[l .8 w h],...
+ 'Callback', @btn_sum_callback);
+uicontrol(commandsp,'Style','pushbutton','String','Clear Outside Sum Channel',...
+ 'Units','normalized', 'Tag', 'btn_clear_outside', ...
+ 'Position',[l .7 w h],...
+ 'Callback', @btn_clear_outside_callback);
+
+uicontrol(commandsp,'Style','pushbutton','String','Detect Spots',...
+ 'Units','normalized', 'Tag', 'btn_detect_spots', ...
+ 'Position',[l .6 w h],...
+ 'Callback', @btn_detect_spots_callback);
+
+uicontrol(commandsp,'Style','pushbutton','String','Cluster Spots',...
+ 'Units','normalized', 'Tag', 'btn_cluster_spots', ...
+ 'Position',[l .5 w h],...
+ 'Callback', @btn_cluster_spots_callback);
+
+uicontrol(commandsp,'Style','pushbutton','String','Reassign Spots',...
+ 'Units','normalized', 'Tag', 'btn_reassign_spots', ...
+ 'Position',[l .4 w h],...
+ 'Callback', @btn_reassign_spots_callback);
+
+% Load WingJ
+uicontrol(commandsp,'Style','pushbutton','String','Load WingJ Data',...
+ 'Units','normalized', 'Tag', 'btn_wingj', ...
+ 'Position',[l .3 w h],...
+ 'Callback', @btn_wingj_callback);
+
+uicontrol(commandsp,'Style','pushbutton','String','Export To CSV',...
+ 'Units','normalized', 'Tag', 'btn_export', ...
+ 'Position',[l .2 w h],...
+ 'Callback', @btn_export_callback);
+
+uicontrol(commandsp,'Style','text',...
+ 'String','Settings',...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'center', ...
+ 'Position',[l .10 w 0.05]);
+
+uicontrol(commandsp,'Style','pushbutton','String','Save',...
+ 'Units','normalized', 'Tag', 'btn_save', ...
+ 'Position',[l .05 w/2 h/2],...
+ 'Callback', @savebtn_callback);
+
+uicontrol(commandsp,'Style','pushbutton','String','load',...
+ 'Units','normalized', 'Tag', 'btn_load', ...
+ 'Position',[l+w/2 .05 w/2 h/2],...
+ 'Callback', @loadbtn_callback);
+
+% Make all the settings
+% Settings Tab
+tgroup = uitabgroup('Parent', settingsp, 'Tag', 'tab_channels');
+
+makeChanSumTab(tgroup);
+makeCleanChanTab(tgroup);
+makeSpotDetectionTab(tgroup);
+makeClusteringTab(tgroup);
+makeReassignTab(tgroup);
+makeWingJTab(tgroup)
+%makeExportTab(tgroup);
+
+end
+
+function makeChanSumTab(parent)
+global X;
+tabTag = 'tab_chan_sum';
+tabName = 'Sum Channels';
+thisTab = findobj('Tag',tabTag);
+if isempty(thisTab)
+ thisTab = uitab('Parent', parent, 'Title', tabName, ...
+ 'Tag', tabTag);
+
+ % Add all the controls
+ %Control settings
+ %Channels of interest as checkboxes
+ sp = 0.05;
+ nChan = X.GetSize('C');
+ uicontrol(thisTab,'Style','text',...
+ 'String','Channels to Sum',...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'left', ...
+ 'Position',[.05 0.9 .3 sp]);
+
+ for k=1:nChan;
+ cbh(k) = uicontrol(thisTab, 'Style','checkbox','String',sprintf('Channel %d',k), ...
+ 'Tag', sprintf('is_sum_chan_%d',k), ...
+ 'Value',1, ...
+ 'Units', 'Normalized', ...
+ 'Position',[.05 .9-sp*k .3 sp]);
+ end
+
+ % Ask if we perform normalization
+ uicontrol(thisTab, 'Style','checkbox','String','Normalize Channels Before Sum', ...
+ 'Tag', 'is_chan_norm', ...
+ 'Value',1, ...
+ 'Units', 'Normalized', ...
+ 'Position',[.05 0.9-sp*(nChan+3) .5 sp]);
+
+end
+end
+
+function makeCleanChanTab(parent)
+tabTag = 'tab_clean';
+tabName = 'Clean Channel';
+thisTab = findobj('Tag',tabTag);
+if isempty(thisTab)
+ thisTab = uitab('Parent', parent, 'Title', tabName, ...
+ 'Tag', tabTag);
+
+ % Add all the controls
+ %Control position
+ sp = 0.05;
+ py = 0.89;
+
+ % Channel for surface creation
+ uicontrol(thisTab,'Style','text',...
+ 'String','Surface Channel',...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py .3 sp]);
+
+ uicontrol(thisTab,'Style','edit',...
+ 'String','',...
+ 'Units','Normalized', ...
+ 'Tag', 'txt_channel_surface', ...
+ 'Position',[.36 py+.01 .1 sp]);
+
+ % Smoothing Factor for Surface
+ uicontrol(thisTab,'Style','text',...
+ 'String','Surface Smoothing',...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py-sp .3 sp]);
+
+ uicontrol(thisTab,'Style','edit',...
+ 'String','5.0',...
+ 'Units','Normalized', ...
+ 'Tag', 'txt_surface_smooth', ...
+ 'Position',[.36 py+.01-sp .1 sp]);
+
+ % Surface Intensity Threshold
+ uicontrol(thisTab,'Style','text',...
+ 'String','Intensity Threshold',...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py-2*sp .3 sp]);
+
+ uicontrol(thisTab,'Style','edit',...
+ 'String','15',...
+ 'Units','Normalized', ...
+ 'Tag', 'txt_surface_intensity_threshold', ...
+ 'Position',[.36 py+.01-2*sp .1 sp]);
+end
+end
+
+function makeSpotDetectionTab(parent)
+tabTag = 'tab_detect';
+tabName = 'Detect Spots';
+thisTab = findobj('Tag',tabTag);
+if isempty(thisTab)
+ thisTab = uitab('Parent', parent, 'Title', tabName, ...
+ 'Tag', tabTag);
+
+ % Add all the controls
+ %Control position
+ sp = 0.05;
+ py = 0.89;
+ % Spot Detection Channel
+ uicontrol(thisTab,'Style','text',...
+ 'String','Spot Channel',...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py .3 sp]);
+
+ uicontrol(thisTab,'Style','edit',...
+ 'String','',...
+ 'Units','Normalized', ...
+ 'Tag', 'txt_channel_spot', ...
+ 'Position',[.36 py+.01 .1 sp]);
+
+ % Spot Spot Sizes
+ uicontrol(thisTab,'Style','text',...
+ 'String','Spot Size [um]',...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py-sp .3 sp]);
+ % XY
+ uicontrol(thisTab,'Style','edit',...
+ 'String','2.8',...
+ 'Units','Normalized', ...
+ 'Tag', 'txt_spot_size_xy', ...
+ 'Position',[.36 py+.01-sp .1 sp]);
+ % Z (Optional)
+ uicontrol(thisTab,'Style','edit',...
+ 'String','2.8',...
+ 'Units','Normalized', ...
+ 'Tag', 'txt_spot_size_z', ...
+ 'Position',[.47 py+.01-sp .1 sp]);
+
+ % Spot Quality Threshold
+ uicontrol(thisTab,'Style','text',...
+ 'String','Spot Quality Threshold',...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py-2*sp .30 sp]);
+
+ uicontrol(thisTab,'Style','edit',...
+ 'String','15',...
+ 'Units','Normalized', ...
+ 'Tag', 'txt_spot_quality_threshold', ...
+ 'Position',[.36 py+.01-2*sp .1 sp]);
+
+ % Remove Spots Closer than
+ uicontrol(thisTab,'Style','checkbox',...
+ 'String','Remove Spots if closer than',...
+ 'Tag', 'is_remove_spots', ...
+ 'Units','Normalized', ...
+ 'Position',[.05 py-3*sp .40 sp]);
+
+ uicontrol(thisTab,'Style','edit',...
+ 'String','1.0',...
+ 'Units','Normalized', ...
+ 'Tag', 'txt_spot_proximity_threshold', ...
+ 'Position',[.45 py-3*sp .07 sp]);
+
+ uicontrol(thisTab,'Style','text',...
+ 'String','um to image edge',...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'left', ...
+ 'Position',[.53 py-3*sp-.01 .30 sp]);
+
+
+end
+end
+
+function makeClusteringTab(parent)
+tabTag = 'tab_cluster';
+tabName = 'Cluster';
+thisTab = findobj('Tag',tabTag);
+if isempty(thisTab)
+ thisTab = uitab('Parent', parent, 'Title', tabName, ...
+ 'Tag', tabTag);
+
+ % Add all the controls
+ % Need 3 tabs: Manual Clustering, Brute Force and KMeans
+ clusMode = uitabgroup('Parent', thisTab, 'Tag', 'tab_cluster_modes');
+
+ %Settings for clustering tab
+ makeSettingsClusterTab(clusMode);
+
+ % Manual Tab
+ makeManualClusterTab(clusMode);
+
+ % Brute Force
+ makeBruteClusterTab(clusMode);
+
+ % KMeans Tab
+ makeKMeansClusterTab(clusMode);
+end
+end
+
+function makeManualClusterTab(parent)
+tabTag = 'tab_cluster_manual';
+tabName = 'Manual';
+thisTab = findobj('Tag',tabTag);
+sp = .05;
+py = .89;
+
+if isempty(thisTab)
+ thisTab = uitab('Parent', parent, 'Title', tabName, ...
+ 'Tag', tabTag);
+
+ uicontrol(thisTab,'Style','checkbox',...
+ 'String','Use Manual Clustering',...
+ 'Tag', 'is_do_manual', ...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py+sp .4 sp]);
+
+
+ uicontrol(thisTab,'Style','text',...
+ 'String','Number of Categories',...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py-sp .3 sp]);
+
+ uicontrol(thisTab,'Style','edit',...
+ 'String','2',...
+ 'Units','Normalized', ...
+ 'Tag', 'txt_manual_cat_num', ...
+ 'Position',[.36 py+.01-sp .1 sp], ...
+ 'Callback', @txt_manual_cat_num_callback);
+
+end
+
+updateManualTab(thisTab);
+end
+
+function txt_manual_cat_num_callback(hObject, eventdata, handles)
+parent = findobj('Tag', 'tab_cluster_manual');
+updateManualTab(parent);
+end
+
+function updateManualTab(parent)
+global X;
+% Add all the controls
+% Grab all the channels that we need to build the settings for
+chanIdx = getChannelIndexes();
+nChan = X.GetSize('C');
+numcat = str2num(get(findobj('Tag', 'txt_manual_cat_num'), 'String'));
+py = 0.89;
+sp = 0.05;
+
+% Make the labels for the columns but only if they do not exist already...
+if isempty(findobj('Tag', 'txt_bg'))
+ uicontrol(parent,'Style','text',...
+ 'String','BG',...
+ 'Units','Normalized', ...
+ 'Tag', sprintf('txt_bg') , ...
+ 'HorizontalAlignment', 'center', ...
+ 'Position',[.26 py-2*sp .1 sp]);
+end
+
+for n=1:(numcat-1)
+ l = findobj('Tag', sprintf('txt_thr_%d', n));
+ if isempty(l)
+ thrtxt(n) = uicontrol(parent,'Style','text',...
+ 'String',sprintf('Thr %d', n),...
+ 'Tag',sprintf('txt_thr_%d', n), ...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'center', ...
+ 'Position',[.26+n*.1 py-2*sp .1 sp], ...
+ 'Visible', 'on');
+ else
+ thrtxt(n) = l;
+ end
+ set(thrtxt(n),'Visible', 'on');
+
+end
+
+
+for k=1:nChan
+ if (any(chanIdx==k))
+ isEnabled = 'on';
+ else
+ isEnabled = 'off';
+ end
+
+ l = findobj('Tag', sprintf('txt_chan_%d', k));
+ if isempty(l)
+ txt(k) = uicontrol(parent,'Style','text',...
+ 'String',sprintf('Channel %d', k),...
+ 'Units','Normalized', ...
+ 'Tag', sprintf('txt_chan_%d',k) , ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py-sp*(k+2) .2 sp]);
+ else
+ txt(k) = l;
+ end
+ set(txt(k),'Enable', isEnabled);
+
+ l= findobj('Tag', sprintf('txt_manual_bg_ch_%d',k));
+ if isempty(l)
+
+ edt(k) = uicontrol(parent,'Style','edit',...
+ 'String','',...
+ 'Units','Normalized', ...
+ 'Tag', sprintf('txt_manual_bg_ch_%d',k) , ...
+ 'Position',[.26 py-sp*(k+2) .1 sp]);
+ else
+ edt(k) = l;
+ end
+ set(edt(k),'Enable', isEnabled);
+
+
+
+ for n=1:(numcat-1)
+ l = findobj('Tag', sprintf('txt_manual_cat_%d_ch_%d',n,k));
+ if isempty(l)
+ thr((k-1)*(numcat-1)+n) = uicontrol(parent,'Style','edit',...
+ 'String','',...
+ 'Units','Normalized', ...
+ 'Tag', sprintf('txt_manual_cat_%d_ch_%d',n,k) , ...
+ 'Position',[.26+n*.1 py-sp*(k+2) .1 sp], ...
+ 'Visible', 'on');
+ else
+ thr((k-1)*(numcat-1)+n)= l;
+ end
+ set(thr((k-1)*(numcat-1)+n),'Enable', isEnabled);
+ set(thr((k-1)*(numcat-1)+n),'Visible', 'on');
+
+ end
+
+ % Delete any columns that are no longer useful
+ for n=(numcat):10
+ obj = findobj('Tag', sprintf('txt_manual_cat_%d_ch_%d',n,k));
+ if ~isempty(obj)
+ delete(obj);
+ end
+ obj = findobj('Tag', sprintf('txt_thr_%d', n));
+ if ~isempty(obj)
+ delete(obj);
+ end
+ end
+
+ % Add button for manual measurement
+ uicontrol(parent,'Style','pushbutton','String','Fill values from spots',...
+ 'Units','normalized', 'Tag', 'btn_calculate_values', ...
+ 'Position',[.2 py-sp*(nChan+6) .4 sp*3], ...
+ 'Callback', @btn_calculate_values_callback);
+
+
+end
+end
+
+function makeBruteClusterTab(parent)
+tabTag = 'tab_cluster_brute_force';
+tabName = 'Brute Force';
+thisTab = findobj('Tag',tabTag);
+sp = .05;
+py = .89;
+
+if isempty(thisTab)
+ thisTab = uitab('Parent', parent, 'Title', tabName, ...
+ 'Tag', tabTag);
+
+ uicontrol(thisTab,'Style','checkbox',...
+ 'String','Use Brute Force Clustering',...
+ 'Tag', 'is_do_brute', ...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py+sp .4 sp]);
+
+ % Number of categories
+ uicontrol(thisTab,'Style','text',...
+ 'String','Number of Categories',...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py-sp .3 sp]);
+
+ uicontrol(thisTab,'Style','edit',...
+ 'String','2',...
+ 'Units','Normalized', ...
+ 'Tag', 'txt_brute_cat_num', ...
+ 'Position',[.36 py+.01-sp .1 sp]);
+end
+end
+
+function makeKMeansClusterTab(parent)
+tabTag = 'tab_cluster_kmeans';
+tabName = 'K-Means';
+thisTab = findobj('Tag',tabTag);
+sp = .05;
+py = .89;
+
+if isempty(thisTab)
+ thisTab = uitab('Parent', parent, 'Title', tabName, ...
+ 'Tag', tabTag);
+
+ uicontrol(thisTab,'Style','checkbox',...
+ 'String','Use K-Means clustering',...
+ 'Tag', 'is_do_kmeans', ...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py+sp .4 sp]);
+
+ uicontrol(thisTab,'Style','text',...
+ 'String','Number of clusters',...
+ 'Tag', 'txt_norm_kmeans', ...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py-sp .3 sp]);
+
+ uicontrol(thisTab,'Style','edit',...
+ 'String','15',...
+ 'Units','Normalized', ...
+ 'Tag', 'txt_kmeans_cluster_num', ...
+ 'Position',[.36 py+.01-sp .1 sp]);
+
+ uicontrol(thisTab,'Style','checkbox',...
+ 'String','Separate K-Means per channel',...
+ 'Tag', 'is_sep_kmeans', ...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py-3*sp .5 sp], ...
+ 'Callback', @is_sep_kmeans_callback);
+
+ uicontrol(thisTab,'Style','text',...
+ 'String','Clusters per channel',...
+ 'Tag', 'txt_sep_kmeans', ...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py-4*sp-.01 .3 sp]);
+
+ uicontrol(thisTab,'Style','edit',...
+ 'String','2',...
+ 'Units','Normalized', ...
+ 'Tag', 'txt_kmeans_sep_cluster_num', ...
+ 'Position',[.36 py-4*sp .1 sp]);
+end
+
+end
+
+function is_sep_kmeans_callback(hObject, eventdata, handles)
+%Get object to enable or not
+sep(1) = findobj('Tag','txt_sep_kmeans');
+sep(2) = findobj('Tag','txt_kmeans_sep_cluster_num');
+
+knorm(1) = findobj('Tag','txt_norm_kmeans');
+knorm(2) = findobj('Tag','txt_kmeans_cluster_num');
+
+
+if get(hObject,'Value')
+ %Disable normal KMeans Values
+ knorm(1).Enable = 'off';
+ knorm(2).Enable = 'off';
+
+ %Enable independent KMeans Values
+ sep(1).Enable = 'on';
+ sep(2).Enable = 'on';
+
+else
+ %Disable normal KMeans Values
+ knorm(1).Enable = 'on';
+ knorm(2).Enable = 'on';
+
+ %Disable independent KMeans Values
+ sep(1).Enable = 'off';
+ sep(2).Enable = 'off';
+
+end
+
+end
+
+function makeSettingsClusterTab(parent)
+tabTag = 'tab_cluster_settings';
+tabName = 'Settings';
+thisTab = findobj('Tag',tabTag);
+sp = .05;
+py = .89;
+
+if isempty(thisTab)
+ thisTab = uitab('Parent', parent, 'Title', tabName, ...
+ 'Tag', tabTag);
+ % Allow empty clusters?
+ uicontrol(thisTab,'Style','checkbox',...
+ 'String','No Empty Custers',...
+ 'Tag', 'is_no_empty_clusters', ...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py+sp .4 sp]);
+
+ % Grouping distance
+ uicontrol(thisTab,'Style','text',...
+ 'HorizontalAlignment', 'right', ...
+ 'String',sprintf('%s\n%s','Groups spots','if closer than [um]'),...
+ 'Units','Normalized', ...
+ 'Position',[.05 py-sp .3 sp+0.02]);
+
+ uicontrol(thisTab,'Style','edit',...
+ 'String','8.5',...
+ 'Units','Normalized', ...
+ 'Tag', 'txt_grouping_distance', ...
+ 'Position',[.36 py+.01-sp .1 sp]);
+
+ % Z Compensation
+ uicontrol(thisTab,'Style','text',...
+ 'String','Z Compensation factor',...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py-2*sp .3 sp]);
+
+ uicontrol(thisTab,'Style','edit',...
+ 'String','3',...
+ 'Units','Normalized', ...
+ 'Tag', 'txt_z_compensation_factor', ...
+ 'Position',[.36 py+.01-2*sp .1 sp]);
+end
+
+end
+
+
+function makeReassignTab(parent)
+tabTag = 'tab_reassign';
+tabName = 'Reassign';
+thisTab = findobj('Tag',tabTag);
+sp = .05;
+h= 0.05;
+py = .89;
+
+if isempty(thisTab)
+ thisTab = uitab('Parent', parent, 'Title', tabName, ...
+ 'Tag', tabTag);
+ % Make one spot per cluster, spot number is cluster size?
+ uicontrol(thisTab,'Style','checkbox',...
+ 'String','Make 1 Spot / Subcluster',...
+ 'Tag', 'is_reassign_single_spot', ...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py+sp .4 h]);
+
+ % Grouping distance
+ uicontrol(thisTab,'Style','text',...
+ 'HorizontalAlignment', 'right', ...
+ 'String','Maximum distance for reassignment [um]',...
+ 'Units','Normalized', ...
+ 'Position',[.05 py-sp .5 h]);
+
+ uicontrol(thisTab,'Style','edit',...
+ 'String','8.5',...
+ 'Units','Normalized', ...
+ 'Tag', 'txt_reassign_max_distance', ...
+ 'Position',[.56 py-sp .1 h]);
+
+ % Z Compensation
+ uicontrol(thisTab,'Style','text',...
+ 'String','Maximum intensity distance [AU]',...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py-2*sp .5 h]);
+
+ uicontrol(thisTab,'Style','edit',...
+ 'String','120',...
+ 'Units','Normalized', ...
+ 'Tag', 'txt_reassign_max_intensity', ...
+ 'Position',[.56 py-2*sp .1 h]);
+
+ % Z Compensation
+ uicontrol(thisTab,'Style','text',...
+ 'String','Reassign if smaller than',...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py-3*sp .5 h]);
+
+ uicontrol(thisTab,'Style','edit',...
+ 'String','10',...
+ 'Units','Normalized', ...
+ 'Tag', 'txt_reassign_cluster_min', ...
+ 'Position',[.56 py-3*sp .1 h]);
+end
+end
+
+function makeWingJTab(parent)
+tabTag = 'tab_wing_j';
+tabName = 'WingJ';
+thisTab = findobj('Tag',tabTag);
+sp = .05;
+py = .89;
+h = .05;
+if isempty(thisTab)
+ thisTab = uitab('Parent', parent, 'Title', tabName, ...
+ 'Tag', tabTag);
+
+ % Make one spot per cluster, spot number is cluster size?
+ uicontrol(thisTab,'Style','checkbox',...
+ 'String','Remove Offset',...
+ 'Tag', 'is_wingj_remove_offset', ...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py .4 h]);
+
+ uicontrol(thisTab,'Style','checkbox',...
+ 'String','Flip Y Coordinate',...
+ 'Tag', 'is_wingj_flip_y', ...
+ 'Units','Normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[.05 py-sp .4 h]);
+
+end
+end
+
+%% Callbacks for all buttons
+function btn_refresh_callback(hObject, eventdata, handles)
+
+end
+
+function btn_sum_callback(hObject, eventdata, handles)
+sumChannels();
+end
+
+function btn_clear_outside_callback(hObject, eventdata, handles)
+clearOutside();
+end
+
+function btn_detect_spots_callback(hObject, eventdata, handles)
+detectSpots();
+end
+
+function btn_cluster_spots_callback(hObject, eventdata, handles)
+clusterSpots();
+end
+
+function btn_reassign_spots_callback(hObject, eventdata, handles)
+reassignSpots();
+end
+
+function btn_export_callback(hObject, eventdata, handles)
+exportToCSV();
+end
+
+function savebtn_callback(hObject, eventdata, handles)
+saveSettings();
+end
+
+function loadbtn_callback(hObject, eventdata, handles)
+loadSettings();
+
+end
+
+function btn_calculate_values_callback(hObject, eventdata, handles)
+
+manualAverageInterface();
+
+end
+
+function btn_wingj_callback(hObject, eventdata, handles)
+loadWingJData();
+end
+
+
+
diff --git a/functions/clearOutside.m b/functions/clearOutside.m
new file mode 100644
index 0000000..1c70cb2
--- /dev/null
+++ b/functions/clearOutside.m
@@ -0,0 +1,42 @@
+function clearOutside()
+
+global X;
+
+ % Get the channel to create the surface from
+ channelCleaning = getFromTag('txt_channel_surface', 'Number');
+
+ % Get the channel we'll use the surface on for cleaning
+ sumChannel = X.GetSize('C');
+
+ % Gaussian smoothing value
+ cleaningSmooth = getFromTag('txt_surface_smooth', 'Number');
+
+ % Surface detection Threshold
+ cleaningThr = getFromTag('txt_surface_intensity_threshold', 'Number');
+
+ % Surface filter, just pick a big one.
+ cleaningFilter = '"Volume" above 200000.0 um^3';
+
+ % Give it a good name.
+ surfaceName = sprintf('Surface For Cleaning, Smooth=%0.2f, Thr=%0.2f', cleaningSmooth, cleaningThr);
+
+
+ % Use EasyXT to detect the surface
+ surface = X.DetectSurfaces(channelCleaning, ...
+ 'Smoothing', cleaningSmooth, ...
+ 'Threshold', cleaningThr, ...
+ 'Filter', cleaningFilter, ...
+ 'Name', surfaceName);
+ % Add it to the Imaris Scene.
+ X.AddToScene(surface);
+
+ % We can now use EasyXT function MakeChannelFromSurfaces to clean the channel
+ X.MakeChannelFromSurfaces(surface, ...
+ 'Add Channel', true, ...
+ 'Mask Channel', sumChannel);
+
+ % Add number of the new channel to the GUI (Channel for Spot Detection)
+ setFromTag('txt_channel_spot', X.GetSize('C'));
+
+end
+
diff --git a/functions/clusterDistances.m b/functions/clusterDistances.m
new file mode 100644
index 0000000..5616767
--- /dev/null
+++ b/functions/clusterDistances.m
@@ -0,0 +1,45 @@
+function clusterDistances( spotData, clusterIndexes, zComp, minDist, parent )
+
+global X;
+
+color = createLUT(max(clusterIndexes), 255);
+for clus = 1:max(clusterIndexes)
+
+ % Get indices of spots belonging to cluster clus
+ clusInd = find(clusterIndexes == clus);
+ nSpots(clus) = numel(clusInd);
+
+
+ spotPositions = spotData.posXYZ(clusInd,:);
+ % Cluster could be empty
+ if ~isempty(spotPositions)
+
+ clusterSpotsFolder = X.CreateGroup(sprintf('Cluster %d', clus));
+ X.AddToScene(clusterSpotsFolder, 'Parent', parent);
+ % compensate Z
+ spotPositions(:,3) = spotPositions(:,3) / zComp;
+
+ % Easy to get subclusters by distance
+ if size(spotPositions, 1) > 1
+ T = clusterdata(spotPositions,'criterion', 'distance', 'distance', 'euclidean', 'cutoff', minDist);
+ else
+ T = 1;
+ end
+
+ % Get the IDs for each cluster and make a spot
+ subclusIds = unique(T);
+ for subClus = subclusIds'
+ idx = find(T==subClus);
+ tempSpot = X.CreateSpots(spotData.posXYZ(clusInd(idx),:), spotData.radii(clusInd(idx),1)', 'Name', sprintf('Subcluster %d', subClus));
+ X.SetColor(tempSpot, color(clus,:));
+ X.AddToScene(tempSpot, 'Parent',clusterSpotsFolder);
+ end
+ ns=length(subclusIds);
+
+ disp(sprintf(' Cluster %d: %d Subclusters: %d Spots', clus, ns, nSpots(clus)));
+ end
+end
+
+
+end
+
diff --git a/functions/clusterSpots.m b/functions/clusterSpots.m
new file mode 100644
index 0000000..2e2e89c
--- /dev/null
+++ b/functions/clusterSpots.m
@@ -0,0 +1,160 @@
+function clusterSpots()
+
+global X;
+
+% Get the kind of clustering we will be doing
+is_do_manual = getFromTag('is_do_manual');
+is_do_brute = getFromTag('is_do_brute');
+is_do_kmeans = getFromTag('is_do_kmeans');
+
+% Grab the Spot Object in the Scene
+nS = X.GetNumberOf('Spots');
+if nS > 1
+ [~, nS] = X.SelectDialog('Spots');
+end
+spots = X.GetObject('Type', 'Spots', 'Number', nS);
+
+% Grab the channels we will be using
+chanIdx = getChannelIndexes();
+
+% Grab the stats we need, namely Intensity Mean
+spotsData = X.GetSelectedStatistics(spots, 'Intensity Mean');
+
+% To cluster, we need an nx2 array with
+% col1: spotID
+% col2: Value of selected statistic
+% Get the indices of the chosen statistic
+
+% Build what we need
+selectedIds = spotsData.ids+1;
+selectedChannels = cell2mat(spotsData.factors(:,2));
+selectedValues = spotsData.values;
+
+
+% for all the Ids and all Channels, will fill all values from the chosen statistic
+for i = 1:length(selectedValues)
+ fullSpotsDataForCluster(selectedIds(i),str2num(selectedChannels(i))) = selectedValues(i);
+end
+
+% Finished data for clustering
+spotsClusterData = fullSpotsDataForCluster(:,chanIdx);
+
+% Clear some memory
+fullSpotsDataForCluster = [];
+
+%Spot Grouping
+spotGrouping = getFromTag('txt_grouping_distance', 'Number');
+
+% Z compensation
+zComp = getFromTag('txt_z_compensation_factor', 'Number');
+
+% Grouping Distance for after clustering
+minDist = getFromTag('txt_grouping_distance', 'Number');
+
+if is_do_manual
+ % Just grab all the channels and all the values and then use the
+ % equivalent of the brute-force to get the results
+ % How many categories per channel
+ nCat = getFromTag('txt_manual_cat_num', 'Number');
+ % Make an array of edges for each channel
+ edges = zeros(numel(chanIdx), nCat+2); % +2 to include 0 and Infinity
+ edges(:,end) = Inf;
+
+ % Everything in double loop, first channels
+ for c=1:numel(chanIdx)
+ %get BG
+ edges(c,2) = getFromTag(sprintf('txt_manual_bg_ch_%d',chanIdx(c)), 'Number');
+ % Now for each category
+ for n=1:(nCat-1)
+ edges(c,2+n) = getFromTag(sprintf('txt_manual_cat_%d_ch_%d',n,chanIdx(c)), 'Number');
+ end
+ end
+
+ [edges, clusterIndexes ] = sortManual(spotsClusterData, edges);
+
+ % create a Parent folder with the Cluster Number
+ manualFolder = X.CreateGroup(sprintf('Manual Clustering %d categories per channel, Z Compensation [%d]', nCat, zComp ));
+ X.AddToScene(manualFolder);
+
+
+ clusterDistances(spotsData, clusterIndexes, zComp, minDist, manualFolder);
+ disp('Manual Clustering DONE');
+ edges
+
+end
+
+if is_do_brute
+ % Here we use a very simple method to cluster them :
+ % Check on each channel for the intensities and categorise them. a
+ % binning into N categories for each channel.
+ nCat = getFromTag('txt_brute_cat_num', 'Number');
+ [centroids, clusterIndexes ] = sortBrute(spotsClusterData, nCat, 160, 10);
+
+ % create a Parent folder with the Cluster Number
+ bruteFolder = XT.CreateGroup(sprintf('Brute Force Clustering %d categories per channel, Z Compensation [%d]', nCat, zComp ));
+ XT.AddToScene(bruteFolder);
+
+ clusterDistances(spotsData, clusterIndexes, zComp, minDist, bruteFolder);
+ disp('Brute Force Clustering DONE');
+ centroids
+end
+
+if is_do_kmeans
+ % Use separate clusterings per channel?
+ is_sep_kmeans = getFromTag('is_sep_kmeans');
+
+ if is_sep_kmeans
+ nClusterPerChannel = getFromTag('txt_kmeans_sep_cluster_num', 'Number');
+ opts = statset('Display','final');
+ str = '';
+ for i=1:numel(chanIdx)
+ [clusterIndexes(:,i), centroids(:,i)] = kmeans(spotsClusterData(:,i),nClusterPerChannel,...
+ 'start','plus',...
+ 'replicates',50,...
+ 'emptyaction','drop',...
+ 'Options',opts);
+
+ str = [str '%d'];
+ end
+
+ [C ia ic] = unique(clusterIndexes, 'rows');
+
+ % for the indexes, make a new array from 1 to to numger of unique
+ % rows
+ newIdx = 1:size(C,1);
+
+ % Get the centroids
+ cent = [];
+ for c=1:numel(chanIdx)
+ cent = [cent centroids(C(:,c),c)];
+ end
+ centroids = cent
+ % create a Parent folder with the Cluster Number
+ sepkMeansFolder = X.CreateGroup(sprintf('Separate K-Means With %d Clusters per Channel, Z Compensation [%d]', nClusterPerChannel, zComp ));
+ X.AddToScene(sepkMeansFolder);
+ clusterDistances(spotsData, newIdx(ic), zComp, minDist, sepkMeansFolder);
+
+ else
+
+ nClusters = getFromTag('txt_kmeans_cluster_num', 'Number');
+
+ opts = statset('Display','final');
+
+ [clusterIndexes, centroids] = kmeans(spotsClusterData,nClusters,...
+ 'start','plus',...
+ 'replicates',50,...
+ 'emptyaction','drop',...
+ 'Options',opts);
+
+ % create a Parent folder with the Cluster Number
+ kMeansFolder = X.CreateGroup(sprintf('K-Means With %d Clusters, Z Compensation [%d]', nClusters, zComp ));
+ X.AddToScene(kMeansFolder);
+
+ clusterDistances(spotsData, clusterIndexes, zComp, minDist, kMeansFolder);
+ disp('K-Means Clustering DONE');
+
+ end
+ centroids
+end
+end
+
diff --git a/functions/createLUT.m b/functions/createLUT.m
new file mode 100644
index 0000000..ccb350d
--- /dev/null
+++ b/functions/createLUT.m
@@ -0,0 +1,4 @@
+function col = createLUT(nClusters, maxInt)
+% to avoid that colorTable get back to red.
+col = hsv(round(nClusters*(1+0.25)))* maxInt;
+end
\ No newline at end of file
diff --git a/functions/detectSpots.m b/functions/detectSpots.m
new file mode 100644
index 0000000..6dda0a0
--- /dev/null
+++ b/functions/detectSpots.m
@@ -0,0 +1,39 @@
+function detectSpots()
+global X;
+
+% Get the channel we'll work on
+channelDetection = getFromTag('txt_channel_spot', 'Number');
+
+% Get the diameter of the spots
+spotDiameterXY = getFromTag('txt_spot_size_xy', 'Number');
+spotDiameterZ = getFromTag('txt_spot_size_z', 'Number');
+
+% Get the quality chosen
+spotQuality = getFromTag('txt_spot_quality_threshold', 'Number');
+
+% Set the spot quality. If it's not set, use Imaris' auto quality.
+if isempty(spotQuality)
+ spotFilter='"Quality" above automatic threshold';
+else
+ spotFilter=sprintf('"Quality" above %d', spotQuality);
+end
+is_remove_spots = getFromTag('is_remove_spots');
+
+% Add a filter to remove spots touching the borders.
+if is_remove_spots
+ spotDist = getFromTag('txt_spot_proximity_threshold', 'Number');
+ spotFilter2 = sprintf('"Distance to Image Border XYZ" above %.3f um',spotDist);
+ spotFilter = [spotFilter ' ' spotFilter2];
+end
+
+% Need to detect spots from the new channel we previously created
+newSpots = X.DetectSpots(channelDetection, ...
+ 'Diameter XY', spotDiameterXY, ...
+ 'Diameter Z', spotDiameterZ, ...
+ 'Spots Filter', spotFilter, ...
+ 'Name', sprintf('Detected Spots from Channel %d', channelDetection) ...
+ );
+
+X.AddToScene(newSpots);
+end
+
diff --git a/functions/exportToCSV.m b/functions/exportToCSV.m
new file mode 100644
index 0000000..7727527
--- /dev/null
+++ b/functions/exportToCSV.m
@@ -0,0 +1,176 @@
+function exportToCSV()
+
+global X;
+
+% Imaris has the coordinates system for the dataset offset from
+% the 'Imaris Origin' (0,0,0), let's make sure we subtract the XMin
+% YMin and ZMin from the Dataset.
+aDataSet = X.ImarisApp.GetDataSet;
+aMinX = aDataSet.GetExtendMinX;
+aMinY = aDataSet.GetExtendMinY;
+aMinZ = aDataSet.GetExtendMinZ;
+
+
+% File for export
+% Which object to export
+groupToExport = X.SelectDialog('Group');
+
+% Returns cell, we just want the string.
+groupToExport= groupToExport{1};
+
+filename = char(X.ImarisApp.GetCurrentFileName);
+[pathstr, name] = fileparts(filename);
+
+% make a save dialog
+[FileName, PathName] = uiputfile('.csv',['Save Results from ' groupToExport],[pathstr name '.csv']);
+
+% Output File
+exportFile = [PathName FileName];
+
+disp('Getting Information on Clusters...');
+nSpots=zeros(0,1);
+subclusterN=zeros(0,1);
+clusterN=zeros(0,1);
+XYZ = zeros(0,3);
+
+%Get info on the currently selected folder...
+folder = X.GetObject('Name', groupToExport);
+
+if strfind(groupToExport, 'Reduced')
+ % We just need the spots and the special statistics that are related to
+ % them.
+ n = X.GetNumberOf('Spots', 'Parent', folder);
+ for i=1:n
+ spots = X.GetObject('Number',i, 'Parent', folder);
+ stat = X.GetSelectedStatistics(spots, 'Number of Spots');
+
+ nS = stat.values;
+ subClN = 1:numel(nS);
+ clN = repmat(i,1,numel(nS));
+ xyz = spots.GetPositionsXYZ;
+
+ nSpots = [nSpots nS'];
+ subclusterN = [subclusterN subClN];
+ clusterN = [clusterN clN];
+ XYZ = [XYZ; xyz];
+ end
+
+else
+
+ %go through children
+ n = X.GetNumberOf('Group', 'Parent', folder);
+
+
+ for i=1:n
+ clus= X.GetObject('Number',i, 'Parent', folder);
+ clusterName = X.GetName(clus);
+ %name: Cluster NUM
+ k = X.GetNumberOf('Spots', 'Parent', clus);
+ for j=1:k
+ spots = X.GetObject('Number',j, 'Parent', clus);
+
+ nSpots(end+1) = size(spots.GetRadii(),1);
+ spotName = X.GetName(spots);
+ subclusterN(end+1) = str2num(spotName(11:end));
+ clusterN(end+1) = str2num(clusterName(8:end));
+
+ xyz = mean(spots.GetPositionsXYZ,1);
+ XYZ(end+1,:) = xyz;
+ end
+
+ end
+end
+
+finalResTable = table(clusterN',subclusterN',nSpots', XYZ, 'VariableNames', {'ClusterNumber', 'SubclusterNumber', 'nSpots', 'XYZ'});
+
+% if WingJ is there, use it... Also use tables
+wingJGrp = X.GetObject('Name', 'WingJ');
+if ~isempty(wingJGrp)
+ nEl = X.GetNumberOf('Spots', 'Parent', wingJGrp);
+
+ markers = X.GetObject('Name', 'Markers');
+ %Get names and positions of these markers
+ mNames = cell(markers.GetNames);
+ mPos = markers.GetPositionsXYZ;
+ mPos = bsxfun(@minus, mPos, [aMinX aMinY aMinZ]);
+
+
+ for i=1:nEl
+ sp = X.GetObject('Type', 'Spots', 'Number', i, 'Parent', wingJGrp);
+ name = X.GetName(sp)
+ posStructure = sp.GetPositionsXYZ;
+ [~, D] = knnsearch(posStructure(:,1:2), XYZ(:,1:2));
+ finalResTable = [finalResTable table(D)];
+
+ k = strfind(name, 'structure');
+ finalResTable.Properties.VariableNames{end} = [strrep(name(k(1):end),'-','_') 'Dist'];
+
+ %Grab the positions for AP
+ if ~isempty(strfind(X.GetName(sp), 'A-P'))
+ posStructure = bsxfun(@minus,sp.GetPositionsXYZ, [aMinX aMinY aMinZ]);
+ SAP = markerSigns(posStructure, XYZ);
+ % Find which marker this corresponds to
+ SmAP = markerSigns(posStructure, mPos);
+
+
+ end
+ %Grab the positions for VD
+ if ~isempty(strfind(X.GetName(sp), 'V-D'))
+ posStructure = bsxfun(@plus,sp.GetPositionsXYZ, [aMinX aMinY aMinZ]);
+ SVD = markerSigns(posStructure, XYZ);
+
+ % Find which marker this corresponds to
+ SmVD = markerSigns(posStructure, mPos);
+ end
+
+ end
+ %Finally Name this thing
+idx1 = find(SAP>0 & SVD>0);
+if ~isempty(idx1)
+ na1 = mNames{find(SmAP>0 & SmVD>0)};
+ locationNames(idx1) = {na1};
+end
+
+idx2 = find(SAP>0 & SVD<0);
+if ~isempty(idx2)
+ na2 = mNames{find(SmAP>0 & SmVD<0)};
+ locationNames(idx2) = {na2};
+end
+
+idx3 = find(SAP<0 & SVD>0);
+if ~isempty(idx3)
+ na3 = mNames{find(SmAP<0 & SmVD>0)};
+ locationNames(idx3) = {na3};
+end
+
+idx4 = find(SAP<0 & SVD<0);
+if ~isempty(idx4)
+ na4 = mNames{find(SmAP<0 & SmVD<0)};
+ locationNames(idx4) = {na4};
+end
+finalResTable.Locations = locationNames';
+
+end
+
+finalResTable.XYZ = bsxfun(@minus, finalResTable.XYZ, [aMinX aMinY aMinZ]);
+
+
+
+
+
+writetable(finalResTable,exportFile);
+
+msgbox(['Results exported to ' exportFile] ,'Export Complete');
+
+% %Prepare table for clustersTable
+% % size(clusterN')
+% % size(subclusterN')
+% % size(nSpots')
+% % size(locationNames')
+% % clusterN'
+% % subclusterN'
+% % nSpots'
+% % locationNames'
+%
+% clusTable = table(clusterN',subclusterN',nSpots', locationNames');
+% clustersTable(groupToExport, table2cell(clusTable));
\ No newline at end of file
diff --git a/functions/getAllData.m b/functions/getAllData.m
new file mode 100644
index 0000000..b5f799c
--- /dev/null
+++ b/functions/getAllData.m
@@ -0,0 +1,89 @@
+function [allIds allSubId allClus allPos allInt allRadii allnSpots] = getAllData()
+% for the current folder that is selected
+% Go through the clusters
+% Get the cluster number
+% Get the subcluster number
+% Get the positions
+% Get the intensities for the needed channels
+
+global X;
+% All data is to be appended
+
+% Which channels for the intensity
+% get from the User Interface the list of channels of interest
+channelsList = getChannelIndexes();
+
+%Get the statistic on which to base clustering
+chosenStat = 'Intensity Mean';
+
+% Z compensation
+zComp = getFromTag('txt_z_compensation_factor', 'Number');
+
+% Which object to Read the data from?
+nFolders = X.GetNumberOf('Group');
+
+if (nFolders > 1)
+ sel = X.SelectDialog('Group');
+else
+ folder = X.GetObject('Type','Group', 'Number', 1);
+ sel = folder.GetName;
+end
+
+mainFolder = X.GetObject('Type','Group', 'Name', sel);
+% Iterate through the clusters
+nClusters = X.GetNumberOf('Group', 'Parent', mainFolder);
+allPos = [];
+allIds = [];
+allInt = [];
+allSub = [];
+allClus = [];
+allnSpots = [];
+allRadii = [];
+
+for i=1:nClusters
+ currCluster = X.GetObject('Type', 'Group', 'Number', i, 'Parent', mainFolder);
+ disp(['Currently Processing Cluster ' num2str(i)]);
+ % Get number of subclusters
+ nSubClus = X.GetNumberOf('Spots', 'Parent', currCluster);
+ clusName = X.GetName(currCluster);
+ clusNum = str2num(clusName(8:end));
+
+ %Iterate through each subcluster
+ for subClus = 1:nSubClus
+
+ %Infer the number of the subcluster, which will give us the subID
+ spot = X.GetObject('Type', 'Spots', 'Number', subClus, 'Parent', currCluster);
+ name = X.GetName(spot);
+ subClusNum = str2num(name(11:end));
+
+ % Get the spots positions
+ posXYZ = spot.GetPositionsXYZ();
+ nSpots = size(posXYZ,1);
+ allPos = [allPos ; posXYZ];
+
+ % Ids from the spots, simple enough?
+ allIds = [allIds 1:nSpots];
+
+ % SubIDs
+ allSub = [allSub; repmat(subClusNum, nSpots,1)];
+
+ % Cluster ID
+ allClus= [allClus; repmat(clusNum, nSpots,1)];
+
+ %Number of spots per subcluster
+ allnSpots = [ allnSpots nSpots ];
+
+ % Now the intensities
+
+ allInt = [allInt; getStatsOnSelectedChannels(spot, chosenStat, channelsList)];
+
+ allRadii = [allRadii; spot.GetRadii()];
+
+
+ end
+end
+
+allSubId = 100000*allClus + allSub;
+
+end
+
diff --git a/functions/getChannelIndexes.m b/functions/getChannelIndexes.m
new file mode 100644
index 0000000..494a5a9
--- /dev/null
+++ b/functions/getChannelIndexes.m
@@ -0,0 +1,10 @@
+function channel_indexes = getChannelIndexes()
+global X;
+nChan = X.GetSize('C');
+channel_indexes =[];
+for c=1:nChan
+ if get(findobj('Tag', sprintf('is_sum_chan_%d',c)),'Value');
+ channel_indexes(end+1) = c;
+ end
+end
+end
\ No newline at end of file
diff --git a/functions/getFromTag.m b/functions/getFromTag.m
new file mode 100644
index 0000000..bccb550
--- /dev/null
+++ b/functions/getFromTag.m
@@ -0,0 +1,20 @@
+function the_thing = getFromTag( tag_name, varargin )
+obj = findobj('Tag', tag_name);
+if isprop(obj,'Style')
+ switch obj.Style
+ case 'edit'
+ the_thing = obj.String;
+ case 'text'
+ the_thing = obj.String;
+ case 'checkbox'
+ the_thing = obj.Value;
+ case 'popupmenu'
+ the_thing = obj.Value;
+
+ end
+ if nargin == 2 && strcmp(varargin{1}, 'Number')
+ the_thing = str2num(the_thing);
+ end
+end
+end
+
diff --git a/functions/getStatsOnSelectedChannels.m b/functions/getStatsOnSelectedChannels.m
new file mode 100644
index 0000000..d696577
--- /dev/null
+++ b/functions/getStatsOnSelectedChannels.m
@@ -0,0 +1,34 @@
+function chStat = getStatsOnSelectedChannels(spotObj, Stat, channelsArray)
+global X;
+spotsData = X.GetSelectedStatistics(spotObj, Stat);
+% Build what we need
+selectedIds = spotsData.ids+1;
+selectedChannels = cell2mat(spotsData.factors(:,2));
+selectedValues = spotsData.values;
+
+for i = 1:length(selectedValues)
+ allStat(selectedIds(i),str2num(selectedChannels(i))) = selectedValues(i);
+end
+
+chStat = allStat(:,channelsArray);
+
+% One spot object per channel in the BG folder
+BGSpotsFolder = X.GetObject('Type', 'Group', 'Name', 'BG');
+if ~isempty(BGSpotsFolder)
+ % The channel number is the name of the spot object
+ nBG = X.GetNumberOf('Spots', 'Parent', BGSpotsFolder);
+ BGVals = zeros(size(channelsArray));
+ for i=1:nBG
+ % Subtract background if it exists
+ BGSpots = X.GetObject('Type', 'Spots', 'Number', i);
+ if ~isempty(BGSpots)
+ BGStats = X.GetSelectedStatistics(BGSpots, 'Intensity Mean');
+ name = X.GetName(BGSpots);
+ channel = str2num(name);
+ disp('Subtracting BG:');
+ BGVals(channel) = BGStats.values(channel,1);
+ end
+ end
+ chStat = bsxfun(@minus,chStat,BGVals);
+end
+
diff --git a/functions/intersections.m b/functions/intersections.m
new file mode 100644
index 0000000..a82ca56
--- /dev/null
+++ b/functions/intersections.m
@@ -0,0 +1,273 @@
+function [x0,y0,iout,jout] = intersections(x1,y1,x2,y2,robust)
+%INTERSECTIONS Intersections of curves.
+% Computes the (x,y) locations where two curves intersect. The curves
+% can be broken with NaNs or have vertical segments.
+%
+% Example:
+% [X0,Y0] = intersections(X1,Y1,X2,Y2,ROBUST);
+%
+% where X1 and Y1 are equal-length vectors of at least two points and
+% represent curve 1. Similarly, X2 and Y2 represent curve 2.
+% X0 and Y0 are column vectors containing the points at which the two
+% curves intersect.
+%
+% ROBUST (optional) set to 1 or true means to use a slight variation of the
+% algorithm that might return duplicates of some intersection points, and
+% then remove those duplicates. The default is true, but since the
+% algorithm is slightly slower you can set it to false if you know that
+% your curves don't intersect at any segment boundaries. Also, the robust
+% version properly handles parallel and overlapping segments.
+%
+% The algorithm can return two additional vectors that indicate which
+% segment pairs contain intersections and where they are:
+%
+% [X0,Y0,I,J] = intersections(X1,Y1,X2,Y2,ROBUST);
+%
+% For each element of the vector I, I(k) = (segment number of (X1,Y1)) +
+% (how far along this segment the intersection is). For example, if I(k) =
+% 45.25 then the intersection lies a quarter of the way between the line
+% segment connecting (X1(45),Y1(45)) and (X1(46),Y1(46)). Similarly for
+% the vector J and the segments in (X2,Y2).
+%
+% You can also get intersections of a curve with itself. Simply pass in
+% only one curve, i.e.,
+%
+% [X0,Y0] = intersections(X1,Y1,ROBUST);
+%
+% where, as before, ROBUST is optional.
+
+% Version: 1.12, 27 January 2010
+% Author: Douglas M. Schwarz
+% Email: dmschwarz=ieee*org, dmschwarz=urgrad*rochester*edu
+% Real_email = regexprep(Email,{'=','*'},{'@','.'})
+
+
+% Theory of operation:
+%
+% Given two line segments, L1 and L2,
+%
+% L1 endpoints: (x1(1),y1(1)) and (x1(2),y1(2))
+% L2 endpoints: (x2(1),y2(1)) and (x2(2),y2(2))
+%
+% we can write four equations with four unknowns and then solve them. The
+% four unknowns are t1, t2, x0 and y0, where (x0,y0) is the intersection of
+% L1 and L2, t1 is the distance from the starting point of L1 to the
+% intersection relative to the length of L1 and t2 is the distance from the
+% starting point of L2 to the intersection relative to the length of L2.
+%
+% So, the four equations are
+%
+% (x1(2) - x1(1))*t1 = x0 - x1(1)
+% (x2(2) - x2(1))*t2 = x0 - x2(1)
+% (y1(2) - y1(1))*t1 = y0 - y1(1)
+% (y2(2) - y2(1))*t2 = y0 - y2(1)
+%
+% Rearranging and writing in matrix form,
+%
+% [x1(2)-x1(1) 0 -1 0; [t1; [-x1(1);
+% 0 x2(2)-x2(1) -1 0; * t2; = -x2(1);
+% y1(2)-y1(1) 0 0 -1; x0; -y1(1);
+% 0 y2(2)-y2(1) 0 -1] y0] -y2(1)]
+%
+% Let's call that A*T = B. We can solve for T with T = A\B.
+%
+% Once we have our solution we just have to look at t1 and t2 to determine
+% whether L1 and L2 intersect. If 0 <= t1 < 1 and 0 <= t2 < 1 then the two
+% line segments cross and we can include (x0,y0) in the output.
+%
+% In principle, we have to perform this computation on every pair of line
+% segments in the input data. This can be quite a large number of pairs so
+% we will reduce it by doing a simple preliminary check to eliminate line
+% segment pairs that could not possibly cross. The check is to look at the
+% smallest enclosing rectangles (with sides parallel to the axes) for each
+% line segment pair and see if they overlap. If they do then we have to
+% compute t1 and t2 (via the A\B computation) to see if the line segments
+% cross, but if they don't then the line segments cannot cross. In a
+% typical application, this technique will eliminate most of the potential
+% line segment pairs.
+
+
+% Input checks.
+error(nargchk(2,5,nargin))
+
+% Adjustments when fewer than five arguments are supplied.
+switch nargin
+ case 2
+ robust = true;
+ x2 = x1;
+ y2 = y1;
+ self_intersect = true;
+ case 3
+ robust = x2;
+ x2 = x1;
+ y2 = y1;
+ self_intersect = true;
+ case 4
+ robust = true;
+ self_intersect = false;
+ case 5
+ self_intersect = false;
+end
+
+% x1 and y1 must be vectors with same number of points (at least 2).
+if sum(size(x1) > 1) ~= 1 || sum(size(y1) > 1) ~= 1 || ...
+ length(x1) ~= length(y1)
+ error('X1 and Y1 must be equal-length vectors of at least 2 points.')
+end
+% x2 and y2 must be vectors with same number of points (at least 2).
+if sum(size(x2) > 1) ~= 1 || sum(size(y2) > 1) ~= 1 || ...
+ length(x2) ~= length(y2)
+ error('X2 and Y2 must be equal-length vectors of at least 2 points.')
+end
+
+
+% Force all inputs to be column vectors.
+x1 = x1(:);
+y1 = y1(:);
+x2 = x2(:);
+y2 = y2(:);
+
+% Compute number of line segments in each curve and some differences we'll
+% need later.
+n1 = length(x1) - 1;
+n2 = length(x2) - 1;
+xy1 = [x1 y1];
+xy2 = [x2 y2];
+dxy1 = diff(xy1);
+dxy2 = diff(xy2);
+
+% Determine the combinations of i and j where the rectangle enclosing the
+% i'th line segment of curve 1 overlaps with the rectangle enclosing the
+% j'th line segment of curve 2.
+[i,j] = find(repmat(min(x1(1:end-1),x1(2:end)),1,n2) <= ...
+ repmat(max(x2(1:end-1),x2(2:end)).',n1,1) & ...
+ repmat(max(x1(1:end-1),x1(2:end)),1,n2) >= ...
+ repmat(min(x2(1:end-1),x2(2:end)).',n1,1) & ...
+ repmat(min(y1(1:end-1),y1(2:end)),1,n2) <= ...
+ repmat(max(y2(1:end-1),y2(2:end)).',n1,1) & ...
+ repmat(max(y1(1:end-1),y1(2:end)),1,n2) >= ...
+ repmat(min(y2(1:end-1),y2(2:end)).',n1,1));
+
+% Force i and j to be column vectors, even when their length is zero, i.e.,
+% we want them to be 0-by-1 instead of 0-by-0.
+i = reshape(i,[],1);
+j = reshape(j,[],1);
+
+% Find segments pairs which have at least one vertex = NaN and remove them.
+% This line is a fast way of finding such segment pairs. We take
+% advantage of the fact that NaNs propagate through calculations, in
+% particular subtraction (in the calculation of dxy1 and dxy2, which we
+% need anyway) and addition.
+% At the same time we can remove redundant combinations of i and j in the
+% case of finding intersections of a line with itself.
+if self_intersect
+ remove = isnan(sum(dxy1(i,:) + dxy2(j,:),2)) | j <= i + 1;
+else
+ remove = isnan(sum(dxy1(i,:) + dxy2(j,:),2));
+end
+i(remove) = [];
+j(remove) = [];
+
+% Initialize matrices. We'll put the T's and B's in matrices and use them
+% one column at a time. AA is a 3-D extension of A where we'll use one
+% plane at a time.
+n = length(i);
+T = zeros(4,n);
+AA = zeros(4,4,n);
+AA([1 2],3,:) = -1;
+AA([3 4],4,:) = -1;
+AA([1 3],1,:) = dxy1(i,:).';
+AA([2 4],2,:) = dxy2(j,:).';
+B = -[x1(i) x2(j) y1(i) y2(j)].';
+
+% Loop through possibilities. Trap singularity warning and then use
+% lastwarn to see if that plane of AA is near singular. Process any such
+% segment pairs to determine if they are colinear (overlap) or merely
+% parallel. That test consists of checking to see if one of the endpoints
+% of the curve 2 segment lies on the curve 1 segment. This is done by
+% checking the cross product
+%
+% (x1(2),y1(2)) - (x1(1),y1(1)) x (x2(2),y2(2)) - (x1(1),y1(1)).
+%
+% If this is close to zero then the segments overlap.
+
+% If the robust option is false then we assume no two segment pairs are
+% parallel and just go ahead and do the computation. If A is ever singular
+% a warning will appear. This is faster and obviously you should use it
+% only when you know you will never have overlapping or parallel segment
+% pairs.
+
+if robust
+ overlap = false(n,1);
+ warning_state = warning('off','MATLAB:singularMatrix');
+ % Use try-catch to guarantee original warning state is restored.
+ try
+ lastwarn('')
+ for k = 1:n
+ T(:,k) = AA(:,:,k)\B(:,k);
+ [unused,last_warn] = lastwarn;
+ lastwarn('')
+ if strcmp(last_warn,'MATLAB:singularMatrix')
+ % Force in_range(k) to be false.
+ T(1,k) = NaN;
+ % Determine if these segments overlap or are just parallel.
+ overlap(k) = rcond([dxy1(i(k),:);xy2(j(k),:) - xy1(i(k),:)]) < eps;
+ end
+ end
+ warning(warning_state)
+ catch err
+ warning(warning_state)
+ rethrow(err)
+ end
+ % Find where t1 and t2 are between 0 and 1 and return the corresponding
+ % x0 and y0 values.
+ in_range = (T(1,:) >= 0 & T(2,:) >= 0 & T(1,:) <= 1 & T(2,:) <= 1).';
+ % For overlapping segment pairs the algorithm will return an
+ % intersection point that is at the center of the overlapping region.
+ if any(overlap)
+ ia = i(overlap);
+ ja = j(overlap);
+ % set x0 and y0 to middle of overlapping region.
+ T(3,overlap) = (max(min(x1(ia),x1(ia+1)),min(x2(ja),x2(ja+1))) + ...
+ min(max(x1(ia),x1(ia+1)),max(x2(ja),x2(ja+1)))).'/2;
+ T(4,overlap) = (max(min(y1(ia),y1(ia+1)),min(y2(ja),y2(ja+1))) + ...
+ min(max(y1(ia),y1(ia+1)),max(y2(ja),y2(ja+1)))).'/2;
+ selected = in_range | overlap;
+ else
+ selected = in_range;
+ end
+ xy0 = T(3:4,selected).';
+
+ % Remove duplicate intersection points.
+ [xy0,index] = unique(xy0,'rows');
+ x0 = xy0(:,1);
+ y0 = xy0(:,2);
+
+ % Compute how far along each line segment the intersections are.
+ if nargout > 2
+ sel_index = find(selected);
+ sel = sel_index(index);
+ iout = i(sel) + T(1,sel).';
+ jout = j(sel) + T(2,sel).';
+ end
+else % non-robust option
+ for k = 1:n
+ [L,U] = lu(AA(:,:,k));
+ T(:,k) = U\(L\B(:,k));
+ end
+
+ % Find where t1 and t2 are between 0 and 1 and return the corresponding
+ % x0 and y0 values.
+ in_range = (T(1,:) >= 0 & T(2,:) >= 0 & T(1,:) < 1 & T(2,:) < 1).';
+ x0 = T(3,in_range).';
+ y0 = T(4,in_range).';
+
+ % Compute how far along each line segment the intersections are.
+ if nargout > 2
+ iout = i(in_range) + T(1,in_range).';
+ jout = j(in_range) + T(2,in_range).';
+ end
+end
+
+% Plot the results (useful for debugging).
+% plot(x1,y1,x2,y2,x0,y0,'ok');
diff --git a/functions/loadSettings.m b/functions/loadSettings.m
new file mode 100644
index 0000000..0269b83
--- /dev/null
+++ b/functions/loadSettings.m
@@ -0,0 +1,32 @@
+function loadSettings()
+[FileName,PathName,FilterIndex] = uigetfile('*.m','Load Settings');
+D = load([PathName FileName], '-mat');
+if ~isfield(D,'data')
+ return
+end
+D.data
+% Restore all you got
+f = findobj('Tag', 'clustetring_main_f');
+%Edit buttons
+txt = findall(f, 'Style', 'edit');
+chk = findall(f, 'Style', 'checkbox');
+rad = findall(f, 'Style', 'radiobutton');
+controls = [txt; chk; rad];
+
+for i=1:numel(D.data)
+ obj = findobj('Tag', D.data(i).tag);
+ if ~isempty(obj)
+ if ~isempty(D.data(i).string)
+ obj.String = D.data(i).string;
+ end
+ if ~isempty(D.data(i).value)
+ obj.Value = D.data(i).value;
+ end
+
+ obj.Visible = D.data(i).visible;
+ end
+
+end
+
+
+end
\ No newline at end of file
diff --git a/functions/loadWingJData.m b/functions/loadWingJData.m
new file mode 100644
index 0000000..87de96f
--- /dev/null
+++ b/functions/loadWingJData.m
@@ -0,0 +1,102 @@
+function loadWingJData()
+
+
+isRemoveOffset = getFromTag('is_wingj_remove_offset');
+isFlipY = getFromTag('is_wingj_flip_y');
+
+% Ask for files to load
+[ filename, pathname ] = uigetfile('*.txt','Select Files to Load', 'MultiSelect', 'on');
+if ~iscell(filename)
+ disp('No files selected for WingJ Import');
+ return
+end
+
+global X;
+aDataSet = X.ImarisApp.GetDataSet;
+[ xy ]= X.GetVoxelSize();
+aMinX = aDataSet.GetExtendMinX;
+aMinY = aDataSet.GetExtendMinY;
+aMaxY = aDataSet.GetExtendMaxY;
+aMinZ = aDataSet.GetExtendMinZ;
+
+% Look for WingJ Group
+parent = X.GetObject('Type', 'Group', 'Name', 'WingJ');
+if ~isempty(parent)
+ X.RemoveFromScene(parent);
+ parent = [];
+end
+
+%Prepare WingJ Folder
+parent = X.CreateGroup('WingJ');
+
+% For each file, use the file name to decide what it is.
+for i=1:numel(filename)
+
+ f = dlmread([pathname filename{i}],'\t');
+ % Ignore Z
+ f(:,3) = 0;
+
+ %Get name of file
+ [path name] = fileparts(filename{i});
+
+ % Calibrate from pixels to microns
+ fCal = f.*xy;
+
+ %Remove offset
+ if isRemoveOffset
+ fCal = bsxfun(@minus,f.*xy,[aMinX aMinY aMinZ]);
+ end
+
+ %Flip Y
+ if isFlipY
+ fCal(:,2) = aMaxY+aMinY - fCal(:,2);
+ end
+
+ if findstr(name, 'A-P')
+ posAP = fCal;
+ end
+
+ if findstr(name, 'V-D')
+ posVD = fCal;
+ end
+
+ spot = X.CreateSpots(fCal, 1, 'Name', name);
+ X.AddToScene(spot, 'Parent', parent);
+
+end
+
+% Now try to find V-D A-P and add these as Measurement Points
+% Find center intersection
+
+[xc, yc, APPos, VDPos] = intersections(posAP(:,1), posAP(:,2), posVD(:,1), posVD(:,2), true );
+xc = xc(1);
+yc = yc(1);
+
+%Find position along curves
+idxCenAP = round(APPos(1));
+idxCenVD = round(VDPos(1));
+%idxCenAP = find(posAP(:,1) == xc & posAP(:,2) == yc)
+%idxCenVD = find(posVD(:,1) == xc & posVD(:,2) == yc)
+
+
+posA = [ posAP(round(idxCenAP/2),1) posAP(round(idxCenAP/2),2) 0];
+posP = [ posAP(round((size(posAP,1)-idxCenAP)/2)+idxCenAP,1) posAP(round((size(posAP,1)-idxCenAP)/2)+idxCenAP,2) 0];
+
+
+posV = [ posVD(round(idxCenVD/2),1) posVD(round(idxCenVD/2),2) 0];
+posD = [ posVD(round((size(posVD,1)-idxCenVD)/2)+idxCenVD,1) posVD(round((size(posVD,1)-idxCenVD)/2)+idxCenVD,2) 0];
+
+% Make the 4 points
+names = {'VA' 'DA' 'VP' 'DP'};
+%test = [posA+posV; posD+posA; posV+posP; posD+posP];
+pos = [posA+posV; posD+posA; posV+posP; posD+posP] - repmat([xc yc 0],4,1);
+
+% Add measurement points
+mp = X.CreateMeasurementPoints(pos, 'Point Names', names, 'Name', 'Markers');
+
+X.AddToScene(mp, 'Parent', parent);
+X.AddToScene(parent);
+
+
+end
+
diff --git a/functions/manualAverageInterface.m b/functions/manualAverageInterface.m
new file mode 100644
index 0000000..e452d6a
--- /dev/null
+++ b/functions/manualAverageInterface.m
@@ -0,0 +1,96 @@
+function manualAverageInterface
+global X;
+global intensities;
+f = figure('Name', 'BIOP - Spots Stats', ...
+ 'NumberTitle', 'off', 'Visible', 'off', 'Position', [0 0 300 200], ...
+ 'Tag', 'main_f', ...
+ 'DockControls', 'off', 'ToolBar', 'none', 'MenuBar', 'none');
+
+% Move the window to the center of the screen.
+movegui(f,'center');
+
+% Make the window visible.
+f.Visible = 'on';
+
+l = .05;
+w = .4;
+h = .1;
+s = 0.1;
+% Refresh button
+uicontrol(f,'Style','pushbutton','String','Refresh',...
+ 'Units','normalized', 'Tag', 'btn_stats_refresh', ...
+ 'Position',[0.10 .75 0.8 h*2],...
+ 'Callback', @btn_stats_refresh_callback);
+% Info text
+uicontrol(f,'Style','text','String','Spots Selected: ',...
+ 'Units','normalized', ...
+ 'HorizontalAlignment', 'right', ...
+ 'Position',[l .60 0.32 h]);
+uicontrol(f,'Style','text','String','',...
+ 'Units','normalized', 'Tag', 'txt_spot_name', ...
+ 'HorizontalAlignment', 'left', ...
+ 'Position',[l+0.33 .60 0.5 h]);
+
+uicontrol(f,'Style','text','String','',...
+ 'Units','normalized', 'Tag', 'txt_spot_stats', ...
+ 'FontWeight', 'bold', ...
+ 'HorizontalAlignment', 'center', ...
+ 'Position',[0.05 .45 0.9 h]);
+
+% 'Add to' popup list
+options = {};
+nChan = getChannelIndexes();
+manCat = getFromTag('txt_manual_cat_num', 'Number');
+for i=nChan
+ options{end+1} = sprintf('Channel %d BG',i);
+ for j=1:(manCat-1)
+ options{end+1} = sprintf('Channel %d Thr %d',i,j);
+ end
+end
+
+
+uicontrol(f, 'Style', 'popup',...
+ 'Units','normalized', 'Tag', 'popup_add_to', ...
+ 'String', options,...
+ 'Position',[l .15 0.5 h]);
+
+% Add to button
+uicontrol(f,'Style','pushbutton','String','Add',...
+ 'Units','normalized', 'Tag', 'btn_add', ...
+ 'Position',[l+0.55 .15 0.2 .1],...
+ 'Callback', @btn_add_callback);
+end
+
+function btn_stats_refresh_callback(hObject, eventdata, handles)
+global X;
+global intensities;
+chan = getFromTag('txt_channel_spot', 'Number');
+spots = X.GetObject();
+idx = spots.GetSelectedIndices();
+stats = X.GetSelectedStatistics(spots, 'Intensity Mean', 'Channel', chan);
+intensities = stats.values(idx+1);
+setFromTag('txt_spot_name', X.GetName(spots));
+setFromTag('txt_spot_stats', sprintf('%d Spots, Average Intensity %.1f', numel(intensities), mean(intensities)));
+
+
+end
+
+function btn_add_callback(hObject, eventdata, handles)
+% Get popup options
+global intensities;
+
+choices = get(findobj('Tag', 'popup_add_to'),'String');
+choice = choices{getFromTag('popup_add_to')};
+
+expression = ['Channel (?<channel>\d+) [BGThr ]+(?<cat>[\d]*).*'];
+tokenNames = regexp(choice,expression,'names');
+
+if isempty(tokenNames.cat)
+ %BG
+ setFromTag(sprintf('txt_manual_bg_ch_%s',tokenNames.channel(1)), sprintf('%.1f',mean(intensities)));
+else
+ setFromTag(sprintf('txt_manual_cat_%s_ch_%s',tokenNames.cat(1), tokenNames.channel(1)), sprintf('%.1f',mean(intensities)));
+end
+
+
+end
\ No newline at end of file
diff --git a/functions/reassignSpots.m b/functions/reassignSpots.m
new file mode 100644
index 0000000..9fcbeea
--- /dev/null
+++ b/functions/reassignSpots.m
@@ -0,0 +1,222 @@
+function [ output_args ] = reassignSpots()
+
+% TODO need to grab all this data from imaris and not from the handles...
+global X;
+
+
+[ ~, allSubId, allClus, allPos, allInt, allRadii, allnSpots ] = getAllData();
+allSubIdNew = allSubId; % New subindexes
+
+maxDist = getFromTag('txt_reassign_max_distance', 'Number');
+maxIntDist = getFromTag('txt_reassign_max_intensity', 'Number');
+isReduced = getFromTag('is_reassign_single_spot');
+
+zComp = getFromTag('txt_z_compensation_factor', 'Number');
+
+% Compensation factor
+intThr = 0.2;
+
+% Grab only the subclusters with 3 spots or less
+uniqueSubIds = unique(allSubId);
+uniqueSubIdsIDX = find(allnSpots <= 3); % subclusters with less than min spot number
+uniqueOtherSubIdsIDX = find(allnSpots > 3); % Other subclusters that are bigger than min spot number
+
+
+% Build an object with only the SPOTS from the uniqueOtherSubIdsIDX indexes
+otherSpotsIDX = [];
+for subIdx = uniqueOtherSubIdsIDX
+ subID = uniqueSubIds(subIdx);
+ %disp(fprintf('SubCluster %d has %d spots', subID, allnSpots(subIdx)));
+ if ~isempty(subID)
+ otherSpotsIDX = [otherSpotsIDX; find(allSubId == subID)];
+ end
+end
+
+% Prepare a 'garbage' cluster so we can add it later.
+garbageSubID = 1;
+garbageID = max(allClus)+1;
+
+% Now work on each unique SubID to get the number of spots, neighbors and
+% reassign them
+for subIdx = uniqueSubIdsIDX
+ % Find the 10 nearest neighbors
+ % First get the position of the spots we are interested in
+ subID = uniqueSubIds(subIdx);
+ if ~isempty(subID)
+ % disp(['Looking at a subcluster of ' num2str(allnSpots(subIdx)) ' spots']);
+
+ % Find the indexes of the spots that belong to this subId
+ spotsIDX = find(allSubId == subID);
+
+ % disp([' Found to be from subclusters ' num2str(allSubId(spotsIDX)')]);
+
+
+ selectedSpots = allPos(spotsIDX,:);
+ candidateSpots = allPos(otherSpotsIDX,:);
+ selectedIntensities = allInt(spotsIDX,:);
+
+ for i=1:size(selectedSpots,1)
+ % For a single spot
+ selectedSpot = selectedSpots(i,:);
+
+ % Z compensation
+ selectedSpotZ = selectedSpot ./ [1 1 zComp];
+ candidateSpotsZ = candidateSpots;
+ candidateSpotsZ(:,3) = candidateSpotsZ(:,3) / zComp;
+
+ % Find k nearest neighbors of a given point
+ % Find the 10 closest points
+ k=10;
+ [IDX,D] = knnsearch(candidateSpotsZ, selectedSpotZ, 'K', k+1);
+ % IDX is mxk
+ % Remove points whose distances are > distThresh
+ usefulIDX = IDX(D < maxDist);
+
+ %Make measurement points to show this
+ %drawMeasurementPoints(XT,candidateSpots(usefulIDX,:),selectedSpotZ )
+
+ %waitforbuttonpress;
+
+ if isempty(usefulIDX)
+ % Do nothing, no reassignment, so put in a
+ % new cluster for 'garbage'
+ allSubIdNew(spotsIDX(i)) = garbageID * 100000 + garbageSubID;
+
+ disp(['Spot with subIDX ' num2str(allSubId(spotsIDX(i))') ' assigned as garbage (Dist)']);
+
+ else
+
+ % Check for valid intensities WRT to this spot
+ % Get the difference of intensities for each channel between the
+ % current spot and all candidates
+
+
+ selectedIntensity = selectedIntensities(i,:);
+ candidateIntensities = allInt(otherSpotsIDX(usefulIDX),:);
+ [IDXInt,DInt] = knnsearch(candidateIntensities, selectedIntensity);
+ % Now apply the intensity tests
+ % Delta of intensity for each channel
+ deltaI = ( candidateIntensities(IDXInt,:) - selectedIntensity).^2;
+ sumDelta = sum(deltaI);
+
+ DIntCorr = DInt .* max(deltaI ./ sumDelta);
+
+ usefulIntIDX = IDXInt((DIntCorr ./intThr) < maxIntDist);
+ if isempty(usefulIntIDX)
+ % Do nothing, no reassignment, so put in a
+ % new cluster for 'garbage'
+ allSubIdNew(spotsIDX(i)) = garbageID * 100000 + garbageSubID;
+ disp(['Spot with subIDX ' num2str(allSubId(spotsIDX(i))') ' assigned as garbage (Int)']);
+
+ else
+ % Ready to reassign spots
+ % Final spot IDs selected from the Distance KNN search
+ subIDFromDist = allSubId(otherSpotsIDX(IDX),:);
+
+ % Final spot IDs selected from the Intensity KNN search
+ subIDFromInt = subIDFromDist(IDXInt,:);
+
+ allSubIdNew(spotsIDX(i)) = subIDFromInt;
+
+
+ disp(['>> Spot with subIDX ' num2str(allSubId(spotsIDX(i))') ' reassigned to subIXD ' num2str(allSubIdNew(spotsIDX(i))')]);
+ disp('>>Distances');
+
+
+ D(D < maxDist)'
+
+ disp('>>Intensities');
+ deltaI
+ DIntCorr
+
+ end
+ end
+ end
+ end
+end
+% finish reassignment using the new subIDs
+
+% Get the new clusters
+clusters = ( allSubIdNew - mod(allSubIdNew,100000) ) / 100000;
+subClus = mod(allSubIdNew,100000);
+
+clusterNumber = max(clusters);
+color = createLUT(clusterNumber, 255);
+
+%Now recreate a group
+if (isReduced)
+ clusterSpotsParentFolder = X.CreateGroup('Reordered Clusters Reduced');
+ X.AddToScene(clusterSpotsParentFolder);
+
+ leftClusters = unique(clusters);
+
+
+ for cId = leftClusters'
+ thePos = zeros(0,3);
+ theRad = zeros(0,3);
+ theNs = zeros(0,1);
+ %Make a spots object for each subcluster average
+ localSubClus = subClus(find(clusters == cId));
+
+ % Get the unique subcluster ID
+ leftSubs = unique(localSubClus);
+
+ for sId = leftSubs'
+ % Find all the spots that are there and create the spot
+ ids = find(clusters == cId & subClus == sId);
+ thePos(end+1,:) = mean(allPos(ids,:),1);
+ theRad(end+1,:) = mean(allRadii((ids),:),1);
+ theNs(end+1) = size(allPos(ids,1),1);
+ end
+ tempSpot = X.CreateSpots(thePos, sqrt(theNs), 'Name', sprintf('Cluster %d', cId));
+ X.AddToScene(tempSpot, 'Parent',clusterSpotsParentFolder);
+ tempColor = color(cId,:);
+ X.SetColor(tempSpot, tempColor);
+
+ % Add the Stat
+ X.AddStatistic(tempSpot, 'Number of Spots', theNs);
+ end
+
+
+
+else
+
+ clusterSpotsParentFolder = X.CreateGroup('Reordered Clusters');
+ X.AddToScene(clusterSpotsParentFolder);
+
+
+ % Find the clusters
+ leftClusters = unique(clusters);
+
+ disp('Reassigning Points in Imaris...');
+
+ % Loop through the clusters
+ for cId = leftClusters'
+ % Make a folder
+ clusterSpotsFolder = X.CreateGroup(sprintf('Cluster %d', cId));
+ X.AddToScene(clusterSpotsFolder, 'Parent', clusterSpotsParentFolder);
+ % Find which subclusters are in this cluster
+ localSubClus = subClus(find(clusters == cId));
+
+ % Get the unique subcluster ID
+ leftSubs = unique(localSubClus);
+
+ % Now create a spot for each subcluster that is here
+ for sId = leftSubs'
+ % Find all the spots that are there and create the spot
+ ids = find(clusters == cId & subClus == sId);
+ s = sprintf('Cluster %d, Subcluster %d has %d spots', cId, sId, size(ids,1));
+ disp(s);
+ tempSpot = X.CreateSpots(allPos((ids),:), allRadii((ids),:), 'Name', sprintf('Subcluster %d', sId));
+ X.AddToScene(tempSpot, 'Parent',clusterSpotsFolder);
+ tempColor = color(cId,:);
+ X.SetColor(tempSpot, tempColor);
+ end
+
+
+ end
+end
+
+disp('Done Reassigning');
+end
+
diff --git a/functions/saveSettings.m b/functions/saveSettings.m
new file mode 100644
index 0000000..c47a034
--- /dev/null
+++ b/functions/saveSettings.m
@@ -0,0 +1,24 @@
+function saveSettings()
+% Controls we need to save are 'edit', 'checkbox', 'radio'
+% so go through all of them and save them in a cell array with the tag and
+% the value or string
+f = findobj('Tag', 'clustetring_main_f');
+%Edit buttons
+txt = findall(f, 'Style', 'edit');
+chk = findall(f, 'Style', 'checkbox');
+rad = findall(f, 'Style', 'radiobutton');
+
+controls = [txt; chk; rad];
+data = struct();
+
+for i=1:numel(controls)
+ data(i).tag = controls(i).Tag;
+ data(i).string = controls(i).String;
+ data(i).value = controls(i).Value;
+ data(i).visible = controls(i).Visible;
+end
+
+% Save this somewhere
+[FileName,PathName,FilterIndex] = uiputfile('settings.m','Save Settings');
+save([PathName FileName],'data');
+end
\ No newline at end of file
diff --git a/functions/setFromTag.m b/functions/setFromTag.m
new file mode 100644
index 0000000..25f1abc
--- /dev/null
+++ b/functions/setFromTag.m
@@ -0,0 +1,14 @@
+function setFromTag( tag_name, value )
+obj = findobj('Tag', tag_name);
+if isprop(obj,'Style')
+ switch obj.Style
+ case 'edit'
+ obj.String = num2str(value);
+ case 'text'
+ obj.String = num2str(value);
+ case 'checkbox'
+ obj.Value = value;
+ end
+ end
+end
+
diff --git a/functions/sortBrute.m b/functions/sortBrute.m
new file mode 100644
index 0000000..bfb7d13
--- /dev/null
+++ b/functions/sortBrute.m
@@ -0,0 +1,24 @@
+function [ centroids, finalIdx ] = sortBrute(data, nCat, theMax, bg)
+
+% Divide the Max by the number of categories
+step = (theMax-bg) / (nCat);
+
+edges = bg:step:(step*nCat);
+edges = [0 edges Inf];
+
+finalIdx = zeros(size(data,1),1);
+
+for i=1:size(data,2)
+ [~, idx] = histc(data(:,i), edges);
+ idx(idx==0) = 1;
+ finalIdx(:,i) = idx;
+
+end
+
+[uClusters] = unique(finalIdx, 'rows');
+
+% Index each row of data.
+
+%Return the centroids
+centroids = (uClusters-1)*step - step/2 + bg;
+centroids(uClusters==1) = bg/2;
\ No newline at end of file
diff --git a/functions/sortManual.m b/functions/sortManual.m
new file mode 100644
index 0000000..b304ba3
--- /dev/null
+++ b/functions/sortManual.m
@@ -0,0 +1,19 @@
+function [ centroids, finalIdx ] = sortManual(data, edges)
+
+for i=1:size(data,2)
+ [~, idx] = histc(data(:,i), edges(i,:));
+ finalIdx(:,i) = idx;
+
+end
+
+[C ia ic] = unique(finalIdx, 'rows');
+
+% for the indexes, make a new array from 1 to to numger of unique
+% rows
+newIdx = 1:size(C,1);
+centroids = [];
+for i=1:size(C,2)
+ centroids = [centroids edges(C(:,i)+2)];
+end
+
+finalIdx = newIdx(ic)';
diff --git a/functions/sumChannels.m b/functions/sumChannels.m
new file mode 100644
index 0000000..9916cef
--- /dev/null
+++ b/functions/sumChannels.m
@@ -0,0 +1,47 @@
+function sumChannels()
+% Get EasyXT
+global X;
+
+is_chan_norm = getFromTag('is_chan_norm');
+
+
+
+% get from the User Interface the list of channels of interest
+channels = getChannelIndexes();
+
+
+% Normalize channels before sum?
+if (is_chan_norm)
+ dataset = X.ImarisApp.GetDataSet();
+
+for i=1:numel(channels)
+ channel = channels(i);
+ chanVol = double(typecast(dataset.GetDataVolumeAs1DArrayBytes(channel-1, 0), 'uint8'));
+ theMin = min(chanVol);
+ theMax = max(chanVol);
+
+ chanVol = uint8(((chanVol-theMin) ./ (theMax-theMin) ) .* 256);
+
+ dataset.SetDataVolumeAs1DArrayBytes(chanVol, channel-1, 0);
+end
+
+
+
+
+% create the string containing channels of interest
+stringToAddChannel = '';
+for index = 1:(numel(channels)-1) %from 1 to ((the length of the Channel's array) -1)
+ chNbr = num2str(channels(index),'%d'); %convert number to string
+ stringToAddChannel = strcat(stringToAddChannel,' ch',chNbr,' + ');% concatenate strings with previous
+end
+
+chNbr = num2str(channels(end),'%d');
+stringToAddChannel = strcat(stringToAddChannel,' ch',chNbr);% string containing channels of interest
+
+% Make sure that it's set as a cell array... Why I am not sure...
+expression = {stringToAddChannel};
+
+sumChannel = X.ChannelMath(expression);
+
+end
+