diff --git a/2_difspot/gtReadAndThreshold.m b/2_difspot/gtReadAndThreshold.m new file mode 100644 index 0000000000000000000000000000000000000000..bf0946ccdfc3a33c3008951fa17f54885925b245 --- /dev/null +++ b/2_difspot/gtReadAndThreshold.m @@ -0,0 +1,99 @@ +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 diff --git a/2_difspot/gtSeedThreshold_doublethr.m b/2_difspot/gtSeedThreshold_doublethr.m index 5382d7b1fbd465ba7065190febda4590eab7cc16..2d5873d619f1d77ef0c79bd9b9da60f664777f69 100644 --- a/2_difspot/gtSeedThreshold_doublethr.m +++ b/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 diff --git a/2_difspot/gtSegmentDiffractionBlobs.m b/2_difspot/gtSegmentDiffractionBlobs.m index 5da50294e9b79426a13fece8bed9327a03448814..90e9ecdd0a0f002aa0074c3cb162c7171a841abf 100644 --- a/2_difspot/gtSegmentDiffractionBlobs.m +++ b/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 diff --git a/2_difspot/gtSegmentDiffractionBlobs_doublethr.m b/2_difspot/gtSegmentDiffractionBlobs_doublethr.m index 49ce98d4f8ee7bd23f92aa7dcc7fe70f16258ac3..e1a6901f63884f99f50d0753ff3af630cd109f55 100644 --- a/2_difspot/gtSegmentDiffractionBlobs_doublethr.m +++ b/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 diff --git a/2_difspot/gtSetupSegmentation.m b/2_difspot/gtSetupSegmentation.m index 7cfaacaeed3a7fbae65794433aeb1c5e52d98d56..45dda432a840f33f168762bce24c03a51e0a248d 100644 --- a/2_difspot/gtSetupSegmentation.m +++ b/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) diff --git a/2_difspot/pfReadAndThreshold.m b/2_difspot/pfReadAndThreshold.m deleted file mode 100644 index edfdb5fe18182c5eaf2a7562f2e30d24289c24ec..0000000000000000000000000000000000000000 --- a/2_difspot/pfReadAndThreshold.m +++ /dev/null @@ -1,100 +0,0 @@ -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 diff --git a/zUtil_DB/dbInsert.m b/zUtil_DB/dbInsert.m index 01491ec73dc1f7de0396633b55565702f4ca59b6..339c36f30358ffd78b014d9ba271b510c1891aa3 100644 --- a/zUtil_DB/dbInsert.m +++ b/zUtil_DB/dbInsert.m @@ -1,20 +1,19 @@ - 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 diff --git a/zUtil_EDF/edf_wait_read.m b/zUtil_EDF/edf_wait_read.m index a1237b313df351bae5e0dcce7ae08a54a87f0305..c76974cf47577fa50753fd3637255b595fec721e 100644 --- a/zUtil_EDF/edf_wait_read.m +++ b/zUtil_EDF/edf_wait_read.m @@ -1,4 +1,4 @@ -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!')