Segmentation/doublethr: speeded up seeds part, by recycling EDF header.

Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>

git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@586 4c865b51-4357-4376-afb4-474e03ccb993

2_difspot/gtReadAndThreshold.m

+function [img, med_value, info] = gtReadAndThreshold(ndx, parameters, thr_single, nodisp, info)
+% GTREADANDTHRESHOLD  Reads the full image with number 'ndx' and thresholds
+% it. It can also return the median value of the full image.
+%
+% It sets the image to zero inside the segmentation bounding box
+% (parameters.seg.bbox) and outside where it is under the threshold.
+% Can subtract the remaining median background of the full image
+% (calculated and applied outside seg.bbox).
+% Thresholds to seg.thr_single by default, or allow another threshold
+% value to be passed in. Pass in NaN to return the greyscale image without
+% thresholding.
+%
+% INPUT
+%   ndx                       - full image number to be read
+%   parameters.
+%     prep.correct_drift
+%     prep.correct_drift_iteration
+%
+%     seg.thr_single          - the threshold value to be applied
+%     seg.bbox                - bounding box of area omitted in segmentation
+%     seg.background_subtract - if true, the full image's median will be
+%                               offset to zero
+%
+% Optional INPUT
+%   thr_single                - threshold value that overrides seg.thr_single;
+%                               if NaN, no thresholding is done
+%   nodisp                    - if true, no text will be displayed
+%   info                      - EDF header structure
+%
+% OUTPUT
+%   im        - segmented image
+%   med_value - median value of the full image
+%   info      - EDF header structure
+%
+% [img, med_value, info] = gtReadAndThreshold(ndx, parameters, thr_single, nodisp, info)
+%
+%
+% Version 002 08-05-2012 by P.Reischig
+%   Updated to new parameters (thr_single) and cleaned.
+%
+% Modified by Nicola Vigano', 2012, nicola.vigano@esrf.fr
+%   Changed gtPlaceSubImage -> gtPlaceSubImage2, added EDF header recycling.
+%
+
+    seg  = parameters.seg;
+    prep = parameters.prep;
+
+    % if a threshold other than th2 is passed in
+    if (~exist('thr_single', 'var'))
+        thr_single = seg.thr_single;
+    end
+    if (~exist('nodisp', 'var'))
+        nodisp = false;
+    end
+    out = GtConditionalOutput(~nodisp);
+
+    % should we wait for images to be shifted by gtApplyDrifts?
+    waitfortranslated = false;
+
+    if isfield(prep, 'correct_drift') && strcmp(prep.correct_drift, 'required')
+        waitfortranslated = prep.correct_drift_iteration;
+    end
+
+    fname = sprintf('1_preprocessing/full/full%04d.edf', ndx);
+
+    out.fprintf('now processing image %s\n', fname);
+
+    if (~exist('info', 'var') || isempty(info))
+        [im_full, info] = edf_wait_read(fname, [], waitfortranslated);
+    else
+        [im_full, info] = edf_wait_read(fname, [], waitfortranslated, info);
+    end
+
+    % would be very easy to add this as an option in gtCreateFullLive - do
+    % it once there and then not again...
+    if (isfield(seg, 'background_subtract') && (seg.background_subtract == true))
+        %subtract the median of the background
+        %taken from gtFullStatCondor - thanks Peter!
+        nullbb = NaN(seg.bbox(4), seg.bbox(3));
+        imn    = gtPlaceSubImage2(nullbb, im_full, seg.bbox(1), seg.bbox(2));
+        imn    = imn(~isnan(imn(:)));
+
+        med_value = median(imn);
+        im_full   = im_full - med_value;
+
+        out.fprintf('subtracting median background value\n');
+    else
+        med_value = NaN;
+    end
+
+    % now blank out the extinction image in the centre
+    img = gtPlaceSubImage2( zeros(seg.bbox(4), seg.bbox(3)), im_full, ...
+                            seg.bbox(1), seg.bbox(2));
+
+    % do the thresholding
+    if ~isnan(thr_single)
+        img(img < thr_single) = 0;
+    end
+end

2_difspot/gtSeedThreshold_doublethr.m

@@ -4,31 +4,29 @@ function gtSeedThreshold_doublethr(first, last, workingdirectory)
 %     --------------------------------------------------------
 %
 %     It finds seeds by thresholding full images individually at the value
-%     parameters.seg.thr_seed. Seeds are the single highest pixels in each 
+%     parameters.seg.thr_seed. Seeds are the single highest pixels in each
 %     connected region.
 %
 %     Puts the seeds into a database table to be used for double threshold
-%     segmention. It also puts the median background value for the full 
+%     segmention. It also puts the median background value for the full
 %     images into a database table.
 %
 %     The segmentation script waits for the full median values table to be
-%     completed before starting. For this reason complete the table even 
+%     completed before starting. For this reason complete the table even
 %     if not subtracting the background, doesn't waste that much effort
 %
 %
+% Version 002 May-2012 by P.Reischig
+%   Update to new parameters names; use of bwconncomp instead of bwlabel for
+%   labelling images; cleaning
 %
+% Modified by Nicola Vigano', 2012, nicola.vigano@esrf.fr
+%   Added EDF header recycling
 %
-%     Version 002 May-2012 by P.Reischig
-%       update to new parameters names; use of bwconncomp instead of bwlabel
-%       for labelling images; cleaning
-%
-
-
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % setup
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 if isdeployed
     first = str2double(first);
     last  = str2double(last);
@@ -41,7 +39,6 @@ end
 fprintf('Changing directory to %s\n',workingdirectory)
 cd(workingdirectory)
 
-
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % load and unpack parameters
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -53,56 +50,52 @@ table_fullmedianvals = sprintf('%sfullmedianvals', parameters.acq.name);
 
 gtDBConnect();
 
-
+info = [];
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % main loop through images
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 for zz = first:last
-    
-    [im, med_val] = pfReadAndThreshold(zz, parameters, NaN);
     % NaN argument means this returns the greyscale image without thresholding
-    
+    [im, med_val, info] = gtReadAndThreshold(zz, parameters, NaN, true, info);
+
     % while doing this, fill in the fullmedianvals table
-    mym(sprintf('update %s set val=%0.20g where ndx=%d', ...
+    mym(sprintf('UPDATE %s SET val=%0.20g WHERE ndx=%d', ...
         table_fullmedianvals, med_val, zz));
-    
+
     % thr_seed threshold
     thr_im = im > parameters.seg.thr_seed;
     thr_im = bwmorph(thr_im, 'open');
-    
+
     % label (2D only)
     % bwconncomp gives:
     %  labels.NumObjects       - number of labels
     %  labels.PixelIdxList{ii} - linear indices in input image for label ii
     labels = bwconncomp(thr_im);
-    
+
     % pick out seeds
     list = [];
-    
+
     for ii = 1:labels.NumObjects
-        
         % marker size test
         if (length(labels.PixelIdxList{ii}) > parameters.seg.seedminarea)
-            
+
             [graylevel, index] = max(im(labels.PixelIdxList{ii}));
-            
+
             [yy, xx] = ind2sub(size(im), labels.PixelIdxList{ii}(index));
-            
+
             list(end+1, 1:5) = [xx yy zz graylevel 0];
             % could also add the bounding boxes of the seeds
         end
     end
-    
-    
+
     % write seeds to database
     % seedID x y z value bad
     if (size(list, 1) > 0)
-        dbInsertMultiple(table_seeds, 'x', list(:,1)', 'y', list(:,2)', ...
-            'z', list(:,3)', 'graylevel', list(:,4)', 'bad', list(:,5)');
+        dbInsertMultiple( table_seeds, ...
+                          'x', list(:,1)', 'y', list(:,2)', 'z', list(:,3)', ...
+                          'graylevel', list(:,4)', 'bad', list(:,5)');
     end
-    
 end % end for
 
-
 end % of main function
 

2_difspot/gtSegmentDiffractionBlobs.m

@@ -6,7 +6,7 @@ function gtSegmentDiffractionBlobs(first, last, workingdirectory)
 % finds and labels connected regions in three dimensions; writes each blob
 % to the database as it finds them to ease database traffic.
 %
-% Uses pfReadAndThreshold that can subtract each full image's remaining 
+% Uses gtReadAndThreshold that can subtract each full image's remaining
 % median background to offset it to zero before segmentation (it is
 % calculated and applied outside the seg.bbox area).
 %
@@ -80,7 +80,7 @@ disp('THIS DOES NOT WRAP OVER THE END OF THE DATASET PROPERLY!')
 
 % read the size of one image to make an estimate for the number we can read
 % at once:
-test    = pfReadAndThreshold(first, parameters, NaN);
+test    = gtReadAndThreshold(first, parameters, NaN);
 details = whos('test');
 
 maxMemoryUsage = 100;  % in MB
@@ -112,7 +112,7 @@ for currentfilendx = first:step:last
     for n = 1:step
         
         % Read and threshold image with parameters.seg.thr_single
-        shortVol(:,:,n) = pfReadAndThreshold(currentfilendx + n - 1, parameters);
+        shortVol(:,:,n) = gtReadAndThreshold(currentfilendx + n - 1, parameters);
         
         % add a median filter to de-noise the thresholded images
         % !!! This a second median filter on the full images that

2_difspot/gtSegmentDiffractionBlobs_doublethr.m

@@ -418,7 +418,7 @@ if all(bb(3:4) > minblobsize(1:2)) % size test
 
         % from original gtSegmentDiffractionBlobs
         mym(dbInsert(table_difblob, 'difblobID', 0))
-        difblobID = mym('select last_insert_id()');
+        difblobID = mym('SELECT last_insert_id()');
         fprintf('***DB WRITE *** Blob %d:inserting %d rows, grown from seed %d\n\n',...
                  difblobID, length(xx), seedID)
 
@@ -565,7 +565,7 @@ for z = bb(3):(bb(3)+bb(6)-1)
     % if we don't have the median value we need, calculate it
     if background_subtract
         if isnan(fullmedianvals(z+1))
-            bg = mym(sprintf('select val from %s where ndx=%d', table_fullmedianvals, z));
+            bg = mym(sprintf('SELECT val FROM %s WHERE ndx=%d', table_fullmedianvals, z));
             fullmedianvals(z+1) = bg;
         else
             bg = fullmedianvals(z+1);
@@ -581,27 +581,12 @@ end
 % set anything inside the seg bounding box to zeros
 % put seg.bbox relative to the volume
 % work out which pixels to set to zero
-firstx = segbb(1)-bb(1)+1;
-firsty = segbb(2)-bb(2)+1;
-lastx  = firstx+segbb(3)-1;
-lasty  = firsty+segbb(4)-1;
-
-% in x
-if firstx < 1
-    firstx = 1;
-end
-if lastx > bb(4)
-    lastx = bb(4);
-end
-
+[firstx, firsty] = deal(segbb(1) - bb(1) + 1, segbb(2) - bb(2) + 1);
+[lastx, lasty] = deal(firstx + segbb(3) - 1, firsty + segbb(4) - 1);
 
-% in y
-if firsty < 1
-    firsty = 1;
-end
-if lasty > bb(5)
-    lasty = bb(5);
-end
+% fix limits
+[firstx, firsty] = max([firstx firsty], 1);
+[lastx, lasty] = min([lastx, lasty], [bb(4), bb(5)]);
 
 if ~(firstx>bb(4) || lastx<1 || firsty>bb(5) || lasty<1)
     % if we have to do something about it

2_difspot/gtSetupSegmentation.m

@@ -331,7 +331,7 @@ sfDoThresholding();
         sfUpdateParameters();
         
         % remove the median background at this point - copy paste from
-        % pfReadAndThreshold
+        % gtReadAndThreshold
         tmpim = im.full(:,:,StackIndex);
         if isfield(parameters.seg, 'background_subtract') && ...
            (parameters.seg.background_subtract==true)

2_difspot/pfReadAndThreshold.m

-function [im, med_value] = pfReadAndThreshold(ndx, parameters, varargin)
-% PFREADANDTHRESHOLD  Reads the full image with number 'ndx' and thresholds
-% it. It can also return the median value of the full image.
-%
-% It sets the image to zero inside the segmentation bounding box
-% (parameters.seg.bbox) and outside where it is under the threshold.
-% Can subtract the remaining median background of the full image 
-% (calculated and applied outside seg.bbox).
-% Thresholds to seg.thr_single by default, or allow another threshold 
-% value to be passed in. Pass in NaN to return the greyscale image without
-% thresholding.
-%
-%
-% INPUT
-%   ndx                       - full image number to be read
-%   parameters.  
-%     prep.correct_drift      
-%     prep.correct_drift_iteration 
-%
-%     seg.thr_single          - the threshold value to be applied
-%     seg.bbox                - bounding box of area omitted in segmentation
-%     seg.background_subtract - if true, the full image's median will be 
-%                               offset to zero
-%
-%   varargin{1}     - threshold value that overrides seg.thr_single;
-%                     if NaN, no thresholding is done
-%   varargin{2}     - if true, no text will be displayed
-%
-%
-% OUTPUT
-%   im        - segmented image
-%   med_value - median value of the full image
-%
-%
-% Version 002 08-05-2012 by P.Reischig
-%   Updated to new parameters (thr_single) and cleaned.
-%
-%
-%
-
-seg  = parameters.seg;
-prep = parameters.prep;
-
-% if a threshold other than th2 is passed in
-if isempty(varargin)
-    thr_single = seg.thr_single;
-else
-    thr_single = varargin{1};
-end
-
-if length(varargin)==2
-    nodisp = varargin{2};
-else
-    nodisp = 0;
-end
-
-% should we wait for images to be shifted by gtApplyDrifts?
-waitfortranslated = false;
-if isfield(prep, 'correct_drift') && strcmp(prep.correct_drift, 'required')
-    waitfortranslated = prep.correct_drift_iteration;
-end
-
-fname = sprintf('1_preprocessing/full/full%04d.edf',ndx);
-
-if ~nodisp
-    disp(sprintf('now processing image %s',fname));
-end
-
-im_full = edf_wait_read(fname, [], waitfortranslated);
-
-% would be very easy to add this as an option in gtCreateFullLive - do
-% it once there and then not again...
-if (isfield(seg, 'background_subtract') && (seg.background_subtract == true))
-    %subtract the median of the background
-    %taken from gtFullStatCondor - thanks Peter!
-    nullbb = NaN(seg.bbox(4), seg.bbox(3));
-    imn    = gtPlaceSubImage(nullbb, im_full, seg.bbox(1), seg.bbox(2));
-    imn    = imn(~isnan(imn(:)));
-    
-    med_value = median(imn);
-    im_full   = im_full - med_value;
-    
-    if ~nodisp
-        disp('subtracting median background value')
-    end
-else
-    med_value = NaN;
-end
-
-% now blank out the extinction image in the centre
-im = gtPlaceSubImage(zeros(seg.bbox(4), seg.bbox(3)), ...
-                     im_full, seg.bbox(1), seg.bbox(2));
-
-                                 
-% do the thresholding
-if ~isnan(thr_single)
-    im(im < thr_single) = 0;
-end
-
-end

zUtil_DB/dbInsert.m

- function cmd=dbInsert(table,varargin)
- % modified 10/12/2007 to use mantissa/exponent (%g) for double precision
- % modified 13/12/2007 to use mantissa/exponent (%0.20g) for double precision /Peter/
+function cmd = dbInsert(table, varargin)
+% modified 10/12/2007 to use mantissa/exponent (%g) for double precision
+% modified 13/12/2007 to use mantissa/exponent (%0.20g) for double precision /Peter/
+% Modified by Nicola Vigano', 2012, ID11 @ ESRF vigano@esrf.eu
 
- beginning=sprintf('insert into %s (',table);
- middle1=[];
- for n=1:2:nargin-1
-   middle1=[middle1 sprintf('%s,',varargin{n})];
- end
- middle1=[middle1(1:end-1) ') values ('];
- middle2=[];
+    titles = varargin(1:2:end);
+    titles(2, :) = {','};
 
- for n=2:2:nargin-1
-   middle2=[middle2 sprintf('"%0.20g",',varargin{n})];
- end
- 
- middle2(end)=')';
- 
- cmd=[beginning middle1 middle2];
+    columns = varargin(2:2:end);
+    numColumns = length(columns);
+    rawColumns = cat(1, columns{:});
 
+    tuplePattern(1, 1:numColumns) = {'%0.20g'};
+    tuplePattern(2, 1:numColumns) = {','};
+    tuplePattern = [tuplePattern{1:end-1}];
+
+    cmd = [ 'INSERT INTO ' table ' (' titles{1:end-1} ') VALUES ' ...
+            '(' sprintf(tuplePattern, rawColumns) ')'];
+end

zUtil_EDF/edf_wait_read.m

-function [img, varargout] = edf_wait_read(fname, bb, waitfortranslated)
+function [img, varargout] = edf_wait_read(fname, bb, waitfortranslated, info)
 % EDF_WAIT_READ.M Reads edf if available, otherwise wait 10s and try again...
 % Usage:
 % data = edf_wait_read(filename[, bb, waitfortranslated]);
@@ -50,7 +50,9 @@ function [img, varargout] = edf_wait_read(fname, bb, waitfortranslated)
     translated = false;
     while ~translated
       try
-        info = edf_info(fname);
+        if (~exist('info', 'var'))
+          info = edf_info(fname);
+        end
         if isfield(info, 'istranslated') && (str2double(info.istranslated) == waitfortranslated)
           translated = true;
           %disp('file has been translated!')