% ss_region_hierarchy( filestem, b_WriteVolume ) % ------------------------------------------------------------------------- % % Component of Scale Space Edges, version 1.2c % % Using over-segmentation results (output of ss_region_average.m), % a binary tree hierarchy is constructed and stored in a Space software % 16bitIndexedTree file format. (The infrastructure for writing each depth % of the tree to an image file is also coded.) % % % % % -> ss_write_Space_16bitIndexedTree.m % % Mark Dow, last modified April 2004 % University of Oregon % Brain Development Lab / Lewis Center for Neuroimaging % http://lcni.uoregon.edu/~mark % dow@uoregon.edu % function ss_region_hierarchy( filestem, b_WriteVolume ) n_VolumePlanesMax = 999999; load( [ filestem '_edge_info' ] ) sz_Image = size( flim_gray_cropped ); % % Generate average gradient of segment image. % flim_ave_gradient = zeros( sz_Image ); % % For every segment: % for i = 1:seg_n % % For every edge point of segment: % for n_ep = seg(i,1) : seg(i,2) % % Mark image at edge point coordinate with average gradient of % % this segment. % if seg(i,5) > flim_ave_gradient( ep(n_ep, 1), ep(n_ep, 2) ) % flim_ave_gradient( ep(n_ep, 1), ep(n_ep, 2) ) = seg(i,5); % % Note: In this function, the average gradient of edges map was only used % % to clean up derived region maps. Marking a binary map with % % each edge point would have been fine. % end % end % end % Make list of sub-region sizes and averages: sub_region(1:n_regions, 1:3 ) = 0; region_coord(1:n_regions ) = 0; region_value(1:n_regions ) = 0; fprintf( '\n Add each labeled region to region list : ' ); % Clean up label map borders for indexed tree index: % Note: flim_regions_all is a copy of the original label map, % with the edge points filled in with neighboring labels. The % original label map, flim_regions_labeled, is modified to % reflect combination of regions. The edge points are not filled % on this map so the number of points within a region is accurately % reflected. But it may be better to assign all points to the nearest % region for value assignments too. flim_regions_all = flim_regions_labeled; flim_regions_all = imdilate( flim_regions_all, strel('diamond',1) ); % Fill blank spots at intersections that dilation missed: % For each pixel: for aa = 1:sz_Image(1) for bb = 1:sz_Image(2) % If blank: if flim_regions_all( aa, bb ) == 0 % If not at right edge: if aa - 1 > 0 % Fill with label to right. flim_regions_all( aa, bb ) = flim_regions_all( aa -1, bb ); % Else: else % Fill with label to left. flim_regions_all( aa, bb ) = flim_regions_all( aa +1, bb ); end end end end for i = 1:n_regions for jj = 1:length( num2str( i - 1 ) ) fprintf( '\b' ); end fprintf( num2str( i ) ); sub_region( i, 1 ) = i; % current super-region ith_region = find( flim_regions_labeled == i ); % coordinates of a region label region_coord(i) = ith_region(1); % a coordinate (for each region) sub_region( i, 2 ) = length( ith_region ); % number in region sub_region( i, 3 ) = flim_average_of_region( ith_region(1) ); % average value of region region_value(i) = flim_average_of_region( ith_region(1) ); % average region value at the 0th depth of tree end % Append segment length to segment list: for i = 1:seg_n seg( i, 11 ) = seg(i,2) - seg(i,1) + 1; end % Make a copy of segment average gradient to use for combined segment average % gradient. seg( : , 12 ) = seg( : , 5 ); % Find regions associated with each segment. % Record a point of each region. fprintf( '\n Find regions across segment: ' ); for i = 1:seg_n for jj = 1:length( num2str( i - 1 ) ) fprintf( '\b' ); end fprintf( num2str( i ) ); % If segment is at least three points long: if seg(i,11) >= 3 seg_nbrs_list( 1:n_regions ) = 0; % For every mid-point (not end-point) for n_ep_mid_seg = seg(i,1) + 1 : seg(i,2) - 1; % If not an edge point if ep(n_ep_mid_seg, 1) - 1 > 0 ... & ep(n_ep_mid_seg, 1) + 1 <= sz_Image(1) ... & ep(n_ep_mid_seg, 2) - 1 > 0 ... & ep(n_ep_mid_seg, 2) + 1 <= sz_Image(2) % Get 3x3 neighborhood of labels. nbr_labels = flim_regions_labeled( ep(n_ep_mid_seg, 1) - 1 : ep(n_ep_mid_seg, 1) + 1, ... ep(n_ep_mid_seg, 2) - 1 : ep(n_ep_mid_seg, 2) + 1 ); % Get list of non-edge point positions (relative to neighborhood), % and labels at those positions. posns = find( nbr_labels > 0 ); nbr_regions = nbr_labels( posns ); % list of neighbor sub-regions, with duplicates % Mark list of sub-regions with a one. for j = 1:length( nbr_regions ) % seg_nbrs_list( nbr_regions(j) ) = 1; % a list with a one at index of region seg_nbrs_list( nbr_regions(j) ) = seg_nbrs_list( nbr_regions(j) ) + 1; % a list with end end end reg_1 = find( seg_nbrs_list == max(seg_nbrs_list) ); n_reg_1 = seg_nbrs_list( reg_1 ); seg_nbrs_list( reg_1(1) ) = 0; reg_2 = find( seg_nbrs_list == max( seg_nbrs_list ) ); n_reg_2 = seg_nbrs_list( reg_2(1) ); % To Do: Find third, fourth region with large number of neighbors % indicating single regions that were pinched off by closely approaching % edges. Figure out which regions should be combined and re-label. if n_reg_1 > 2 & n_reg_2 > 2 & reg_1(1) ~= reg_2(1) seg( i, 9 ) = reg_1(1); seg( i, 10 ) = reg_2(1); else seg(i,12) = 999999; end else seg(i,12) = 999999; end end % Remove segments that have no (=0) average gradient. % To Do: Do any? Why? for k = 1:seg_n if seg(k,5) < .00000001 seg(k,5) = 99999; end if seg(k,12) < .00000001 seg(k,12) = 99999; end end % To Do: One and only one method must be selected. b_AveSegGradSimMet = 0; b_RegValDifSimMet = 0; b_SegValBlendDifSimMet = 1; flim_average_of_region_cmb = flim_average_of_region; v_max_region_value = max(max( flim_average_of_region )); flim_average_of_region_cmb_clean = 0; i_VolumePlane = 0; vol_region_hierarchy( 1:sz_Image(2), 1:sz_Image(1), 1 ) = 0; reg_pairs_info = 0; n_reg_pairs = 1; n_reg_pairs_l = 1; n_segs_of_pairs = 1; n_reg_pairs_col = 5; LUT( 1:65536, 1:256 ) = 0; LUT = uint8(LUT); LUT_column(1:n_regions) = 0; % Get extrema and normalize of full over-segmentation LUT. LUT_column(1:n_regions) = region_value( 1:n_regions ); region_value_min = min(min(LUT_column)); LUT_column = LUT_column - region_value_min; region_value_max = max(max(LUT_column)); LUT_column = floor( 255*(LUT_column/region_value_max) ); % Fill second LUT column with full over-segmentation LUT. First column is % 0, reserved for the empty set. LUT( 2:n_regions+1, 1 ) = uint8(LUT_column'); v_max = max( sub_region( :, 3 ) ); v_min = min( sub_region( :, 3 ) ); ss_max = max( seg( :, 5 ) ); ss_min = min( seg( :, 5 ) ); % While there is still a segment on the list to be eliminated: while length( find( seg( : , 12 ) < 99999 ) ) > 0 & n_reg_pairs > 0 % Find similarity metric between region pairs of each segment: seg_pairs_checked( 1:seg_n ) = 0; reg_pairs_info(1:n_reg_pairs,1:n_reg_pairs_col) = 0; i_rpi = 0; n_reg_pairs = 0; % For each segment: fprintf( '\n Add to region pair list: ' ); for i_spc = 1:seg_n for jj = 1:length( num2str( i_spc - 1 ) ) fprintf( '\b' ); end fprintf( num2str( i_spc ) ); % If segment was removed, don't check its region pair. if seg( i_spc, 12 ) >= 99999 seg_pairs_checked( i_spc ) = 1; end % If pair might be new: if seg_pairs_checked( i_spc ) == 0 % Get the pair. reg_pair(1) = seg( i_spc, 9 ); reg_pair(2) = seg( i_spc, 10 ); seg_pairs_checked( i_spc ) = 1; % If it's a valid pair: if reg_pair(1) > 0 & reg_pair(2) > 0 ... & reg_pair(1) ~= reg_pair(2) i_rpi = i_rpi + 1; n_reg_pairs_l = max( n_reg_pairs_l, i_rpi ); n_reg_pairs = i_rpi; reg_pairs_info( i_rpi, 3 ) = reg_pair(1); % label of first of pair reg_pairs_info( i_rpi, 4 ) = reg_pair(2); % label of second of pair % Pair similarity metric: if b_AveSegGradSimMet == 1 % - average gradient of separation segments reg_pairs_info( i_rpi, 1 ) = seg( i_spc, 5 ); reg_pairs_info( i_rpi, 5 ) = reg_pairs_info( i_rpi, 1 ); end if b_RegValDifSimMet == 1 % - average value difference of regions reg_pairs_info( i_rpi, 1 ) = abs( sub_region( seg( i_spc, 9 ), 3 ) - sub_region( seg( i_spc, 10 ), 3 ) ); reg_pairs_info( i_rpi, 5 ) = reg_pairs_info( i_rpi, 1 ); end if b_SegValBlendDifSimMet == 1 % - average value difference of regions reg_pairs_info( i_rpi, 1 ) = seg( i_spc, 5 ); reg_pairs_info( i_rpi, 5 ) = .5*reg_pairs_info( i_rpi, 1 )/ss_max - ss_min + .5*( abs( sub_region( seg( i_spc, 9 ), 3 ) - sub_region( seg( i_spc, 10 ), 3 ) ) )/( v_max - v_min ); end % Pair separation segments length. reg_pairs_info( i_rpi, 2 ) = seg( i_spc, 11 ); % Check the rest of the list for other segments separating this pair. cc = reg_pair(1); dd = reg_pair(2); for i_spc_after = i_spc + 1 : seg_n % If this segment is a valid segment: if seg( i_spc_after, 12 ) < 99999 % And has the same region pair: aa = seg( i_spc_after, 9 ); bb = seg( i_spc_after, 10 ); if ( ( aa == cc & bb == dd ) | ( aa == dd & bb == cc ) ) seg_pairs_checked( i_spc_after ) = 1; n_sp = reg_pairs_info( i_rpi, 2 ); n_after = seg( i_spc_after, 11 ); % Update pair similarity metric: if b_AveSegGradSimMet == 1 % - average gradient of separation segments reg_pairs_info( i_rpi, 1 ) = ( reg_pairs_info( i_rpi, 1 )*n_sp ... + seg( i_spc_after, 5 )*n_after ) ... / ( n_sp + n_after ); reg_pairs_info( i_rpi, 5 ) = reg_pairs_info( i_rpi, 1 ); end % Update pair similarity metric: if b_RegValDifSimMet == 1 % doesn't change end % Update pair similarity metric: if b_SegValBlendDifSimMet == 1 % - average gradient of separation segments reg_pairs_info( i_rpi, 1 ) = ( reg_pairs_info( i_rpi, 1 )*n_sp ... + seg( i_spc_after, 5 )*n_after ) ... / ( n_sp + n_after ); reg_pairs_info( i_rpi, 5 ) = .5*reg_pairs_info( i_rpi, 1 )/ss_max ... + .5*( abs( sub_region( seg( i_spc, 9 ), 3 ) - sub_region( seg( i_spc, 10 ), 3 ) ) )/( v_max - v_min ); end % Update pair separation segments length. reg_pairs_info( i_rpi, 2) = n_sp + n_after; end end end end end % Take this segment off the list. seg( i_spc, 12 ) = 999999; end end combination_difference = 0; combination_difference_nth_plane = 0; n_combinations = 0; regions_combined_this_plane = 0; i_regions_combined_this_plane = 0; last_n_reg_pairs = n_reg_pairs; fprintf( '\n Region pairs left: ' ); % While there is still a segment on the list to be eliminated: while n_reg_pairs > 0 % Backspace over last number printed. for jj = 1:length( num2str( last_n_reg_pairs ) ) fprintf( '\b' ); end % Print this last_n_reg_pairs = n_reg_pairs; fprintf( num2str( n_reg_pairs ) ); if n_reg_pairs > -1 % Find region pair(s) with smallest similarity metric. i_rp = find( reg_pairs_info( 1:n_reg_pairs, 5 ) == min( reg_pairs_info( 1:n_reg_pairs, 5 ) ) ); n_combinations = n_combinations + 1; combination_difference( n_combinations ) = reg_pairs_info( i_rp(1), 5 ); combination_difference_nth_plane( n_combinations ) = i_VolumePlane; if n_combinations > 1 combination_difference( n_combinations ) = max( combination_difference( n_combinations - 1 ), combination_difference( n_combinations ) ); end % Combine the regions across this segment: % Get number and value of two regions. i_reg_1 = reg_pairs_info( i_rp(1), 3 ); i_reg_2 = reg_pairs_info( i_rp(1), 4 ); n_Reg_1 = sub_region( i_reg_1, 2 ); n_Reg_2 = sub_region( i_reg_2, 2 ); v_Reg_1 = sub_region( i_reg_1, 3 ); v_Reg_2 = sub_region( i_reg_2, 3 ); % If one of the regions is the result of a combination since the % last plane was written: if length( find( regions_combined_this_plane == i_reg_1 ) ) > 0 ... | length( find( regions_combined_this_plane == i_reg_2 ) ) > 0 regions_combined_this_plane = 0; i_regions_combined_this_plane = 0; if b_WriteVolume == 1 ... & ( ( n_reg_pairs < n_VolumePlanesMax ... & i_VolumePlane < n_VolumePlanesMax ) ... | n_reg_pairs <= 1 ) i_VolumePlane = i_VolumePlane + 1; flim_average_of_region_cmb_clean = flim_average_of_region_cmb; % flim_average_of_region_cmb_clean( find( flim_ave_gradient ) ) = 0; % flim_average_of_region_cmb_clean = imdilate( flim_average_of_region_cmb_clean, strel( 'diamond', 1 ) ); % % Fill blank spots at intersections: % % For each pixel: % for aa = 1:sz_Image(1) % for bb = 1:sz_Image(2) % % % If blank: % if flim_average_of_region_cmb_clean( aa, bb ) == 0 % % % If not at right edge: % if aa - 1 > 0 % % % Fill with pixel to right. % flim_average_of_region_cmb_clean( aa, bb ) = flim_average_of_region_cmb_clean( aa -1, bb ); % % % Else: % else % % % Fill with pixel to left. % flim_average_of_region_cmb_clean( aa, bb ) = flim_average_of_region_cmb_clean( aa +1, bb ); % end % end % end % end flim_average_of_region_cmb_clean = flim_average_of_region_cmb_clean / v_max_region_value; % imwrite( flim_average_of_region_cmb_clean, [ filestem '_' num2str(i_VolumePlane) '.tif'] ); region_value( 1:n_regions ) = flim_average_of_region_cmb( region_coord(1:n_regions) ); % average region value at the 0th level of tree % Get extrema and normalize of full over-segmentation LUT. LUT_column(1:n_regions) = region_value( 1:n_regions ); LUT_column = LUT_column - region_value_min; LUT_column = floor( 255*(LUT_column/region_value_max) ); LUT( 2:n_regions+1, i_VolumePlane ) = uint8(LUT_column'); end end i_reg_nbr_pairs_1 = 0; i_reg_nbr_pairs_2 = 0; i_reg_nbr_pairs_3 = 0; i_reg_nbr_pairs_4 = 0; i_reg_nbr_pairs_1 = find( reg_pairs_info(1:n_reg_pairs, 3) == i_reg_1 ); i_reg_nbr_pairs_2 = find( reg_pairs_info(1:n_reg_pairs, 4) == i_reg_1 ); % Append indices for second of pair is first region. i_reg_nbr_pairs_1( length(i_reg_nbr_pairs_1)+1:length(i_reg_nbr_pairs_1)+length(i_reg_nbr_pairs_2) ) = i_reg_nbr_pairs_2; i_reg_nbr_pairs_3 = find( reg_pairs_info(1:n_reg_pairs, 3) == i_reg_2 ); i_reg_nbr_pairs_4 = find( reg_pairs_info(1:n_reg_pairs, 4) == i_reg_2 ); i_reg_nbr_pairs_3( length(i_reg_nbr_pairs_3)+1:length(i_reg_nbr_pairs_3)+length(i_reg_nbr_pairs_4) ) = i_reg_nbr_pairs_4; % If regions of pair are different (they shouldn't ever be the % same, but for any abberant dangling segments): if i_reg_1 ~= i_reg_2 ... & ( 10*reg_pairs_info( i_rp(1), 2 )*reg_pairs_info( i_rp(1), 2 ) > n_Reg_1 ... | 10*reg_pairs_info( i_rp(1), 2 )*reg_pairs_info( i_rp(1), 2 ) > n_Reg_2 ) % Island region preservation condition. % & ( length(i_reg_nbr_pairs_1) > 1 & length(i_reg_nbr_pairs_3) > 1 ) ... % Add the resulting region to list of touched regions since % last plane written. i_regions_combined_this_plane = i_regions_combined_this_plane + 1; regions_combined_this_plane( i_regions_combined_this_plane ) = i_reg_1; total_number_in_two_regions = n_Reg_1 + n_Reg_2; ave_value_of_two_regions = ( n_Reg_1*v_Reg_1 + n_Reg_2*v_Reg_2 ) / total_number_in_two_regions; % Redraw both regions with the average value. flim_average_of_region_cmb( find( flim_regions_labeled == i_reg_1 ) ) = ave_value_of_two_regions; flim_average_of_region_cmb( find( flim_regions_labeled == i_reg_2 ) ) = ave_value_of_two_regions; % Relabel all of second region with first region's label. flim_regions_labeled( find( flim_regions_labeled == i_reg_2 ) ) = i_reg_1; % Update combined region stats. sub_region( i_reg_1, 2 ) = total_number_in_two_regions; sub_region( i_reg_1, 3 ) = ave_value_of_two_regions; % Update boundary info of regions that are adjacent to both % regions to be combined: for nth_pair = 1 : length( i_reg_nbr_pairs_1 ) i_pair_1 = i_reg_nbr_pairs_1(nth_pair); nbr_1 = 0; nbr_2 = 0; % Get index of neighbor of first region. if reg_pairs_info( i_pair_1, 3 ) == i_reg_1 if reg_pairs_info( i_pair_1, 4 ) ~= i_reg_1 nbr_1 = reg_pairs_info( i_pair_1, 4 ); end else if reg_pairs_info( i_pair_1, 3 ) ~= i_reg_1 nbr_1 = reg_pairs_info( i_pair_1, 3 ); end end for nth_pair_2 = 1 : length( i_reg_nbr_pairs_3 ) i_pair_2 = i_reg_nbr_pairs_3(nth_pair_2); % Get index of neighbor of second region. if reg_pairs_info( i_pair_2, 3 ) == i_reg_2 if reg_pairs_info( i_pair_2, 4 ) ~= i_reg_2 nbr_2 = reg_pairs_info( i_pair_2, 4 ); end else if reg_pairs_info( i_pair_2, 3 ) ~= i_reg_2 nbr_2 = reg_pairs_info( i_pair_2, 3 ); end end % If a pair is adjacent to both regions that are combined: if nbr_1 == nbr_2 & nbr_1 ~= 0 & nbr_2 ~= 0 reg_pairs_info( i_pair_1, 1 ) = ... ( reg_pairs_info( i_pair_1, 2 )*reg_pairs_info( i_pair_1, 1 ) ... + reg_pairs_info( i_pair_2, 2 )*reg_pairs_info( i_pair_2, 1 ) ) ... / ( reg_pairs_info( i_pair_1, 2 ) + reg_pairs_info( i_pair_2, 2 ) ); reg_pairs_info( i_pair_1, 2 ) = ... reg_pairs_info( i_pair_1, 2 ) + reg_pairs_info( i_pair_2, 2 ); reg_pairs_info( i_pair_1, 5 ) = ... .5*( reg_pairs_info( i_pair_1, 1 )/ss_max ) ... + .5*( abs( sub_region( reg_pairs_info( i_pair_1, 3 ), 3 ) - sub_region( reg_pairs_info( i_pair_1, 4 ), 3 ) )/( v_max - v_min ) ); reg_pairs_info( i_pair_2, : ) = reg_pairs_info( n_reg_pairs, : ); n_reg_pairs = n_reg_pairs - 1; end end end i_replace_1 = find( reg_pairs_info(1:n_reg_pairs, 3) == i_reg_2 ); for ir = 1:length(i_replace_1) if reg_pairs_info( i_replace_1(ir), 4 ) ~= i_reg_1 reg_pairs_info( i_replace_1(ir), 3 ) = i_reg_1; end end i_replace_2 = find( reg_pairs_info(1:n_reg_pairs, 4) == i_reg_2 ); for ir = 1:length(i_replace_2) if reg_pairs_info( i_replace_2(ir), 3 ) ~= i_reg_1 reg_pairs_info( i_replace_2(ir), 4 ) = i_reg_1; end end end if ~( 10*reg_pairs_info( i_rp(1), 2 )*reg_pairs_info( i_rp(1), 2 ) > n_Reg_1 ... | 10*reg_pairs_info( i_rp(1), 2 )*reg_pairs_info( i_rp(1), 2 ) > n_Reg_2 ) reg_pairs_info( i_rp(1), 1 ) = min( 1.2*reg_pairs_info( i_rp(1), 1 ), ss_max ); end reg_pairs_info( i_rp(1), : ) = reg_pairs_info( n_reg_pairs, : ); n_reg_pairs = n_reg_pairs - 1; end % n_reg_pairs > 0 end % figure % plot(combination_difference); % figure % plot(combination_difference_nth_plane); for i = -7:7 DG_filter(i + 8) = -dgauss(i,2); end % figure % plot(DG_filter); G1_combo = conv( combination_difference, DG_filter ); G1_combo = G1_combo(7:length( G1_combo ) - 7); % figure % plot(G1_combo); % baseline_slope = mean( G1_combo( min(1, round( length(G1_combo)/20) ) : round( length(G1_combo)/2 ) ) ); % b_written = 0; % reference_plane = 0; % % for i = round( length(G1_combo)/2 ) : length(G1_combo) % % if G1_combo(i) > 20*baseline_slope & b_written == 0; % % flim_average_of_region_cmb_clean = imread( [ filestem '_' num2str(combination_difference_nth_plane(i)) '.tif'] ); % % imwrite( flim_average_of_region_cmb_clean, [ filestem '_reg_' num2str( 7 ) '.tif'] ); % vol_region_hierarchy( :, :, 7 ) = flim_average_of_region_cmb_clean'; % % b_written = 1; % reference_plane = combination_difference_nth_plane(i); % end % end % % if reference_plane == 0 % % flim_average_of_region_cmb_clean = imread( [ filestem '_' num2str( i_VolumePlane - 2 ) '.tif'] ); % % imwrite( flim_average_of_region_cmb_clean, [ filestem '_reg_' num2str( 7 ) '.tif'] ); % vol_region_hierarchy( :, :, 7 ) = flim_average_of_region_cmb_clean'; % reference_plane = i_VolumePlane - 2; % end % % for i = 0:5 % n = round( i*reference_plane/6 ) + 1; % flim_average_of_region_cmb_clean = imread( [ filestem '_' num2str( n ) '.tif'] ); % % imwrite( flim_average_of_region_cmb_clean, [ filestem '_reg_' num2str( i + 1 ) '.tif'] ); % vol_region_hierarchy( :, :, i + 1 ) = flim_average_of_region_cmb_clean'; % end % n = round( (i_VolumePlane + reference_plane )/2 ); % flim_average_of_region_cmb_clean = imread( [ filestem '_' num2str( n ) '.tif'] ); % % imwrite( flim_average_of_region_cmb_clean, [ filestem '_reg_' num2str( 8 ) '.tif'] ); % vol_region_hierarchy( :, :, 8 ) = flim_average_of_region_cmb_clean'; % If the tree depth is larger than a byte: if i_VolumePlane > 256 % Use the top of the tree, keeping the full over-segmentation as the leaves of the tree. LUT( 1:n_regions+1, 2:256 ) = LUT( 1:n_regions+1, i_VolumePlane-254:i_VolumePlane ); LUT = LUT(:, 1:256 ); end if b_WriteVolume == 1 ss_write_Space_16bitIndexedTree( flim_regions_all', LUT, [ filestem '_region_tree' ] ); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % y = gauss( x, std ) % ------------------------ function y = gauss( x, std ) y = exp(-x^2/(2*std^2)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % y = dgauss( x, std ) % -------------------------- function y = dgauss( x, std ) y = x * gauss(x,std);