Ignore:
Timestamp:
Jun 19, 2016, 8:32:16 PM (8 years ago)
Author:
sommeria
Message:

ima2vol mofidfied for vertical cuts

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/series/ima2vol.m

    r924 r953  
    1 %'ima2vol': concatene  image series to form a 'volume' image .vol used for 3D PIV
     1%'ima2vol': concatene  image series to form a 'volume' image, make vertical cuts along x and y
    22%------------------------------------------------------------------------
    33% function GUI_input=ima2vol(num_i1,num_i2,num_j1,num_j2,Series)
     
    3333%=======================================================================
    3434
    35 function GUI_input=ima2vol(num_i1,num_i2,num_j1,num_j2,Series)
    36 %requests for the visibility of input windows in the GUI series  (activated directly by the selection in the menu ACTION)
    37 if ~exist('num_i1','var')
    38     GUI_input={};
    39     return %exit the function
     35function ParamOut=ima2vol(Param)
     36
     37%% set the input elements needed on the GUI series when the function is selected in the menu ActionName or InputTable refreshed
     38if isstruct(Param) && isequal(Param.Action.RUN,0)
     39    ParamOut.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
     40    ParamOut.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
     41    ParamOut.NbSlice='off'; %nbre of slices ('off' by default)
     42    ParamOut.VelType='one';% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
     43    ParamOut.FieldName='one';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
     44    ParamOut.FieldTransform = 'on';%can use a transform function
     45    ParamOut.TransformPath=fullfile(fileparts(which('uvmat')),'transform_field');% path to transform functions (needed for compilation only)
     46    ParamOut.ProjObject='on';%can use projection object(option 'off'/'on',
     47    ParamOut.Mask='on';%can use mask option   (option 'off'/'on', 'off' by default)
     48    ParamOut.OutputDirExt='.vertical_cut';%set the output dir extension
     49    ParamOut.OutputFileMode='NbInput';% '=NbInput': 1 output file per input file index, '=NbInput_i': 1 file per input file index i, '=NbSlice': 1 file per slice
     50      %check the input files
     51    ParamOut.CheckOverwriteVisible='on'; % manage the overwrite of existing files (default=1)
     52    first_j=[];
     53    if isfield(Param.IndexRange,'first_j'); first_j=Param.IndexRange.first_j; end
     54    PairString='';
     55    if isfield(Param.IndexRange,'PairString'); PairString=Param.IndexRange.PairString; end
     56    [i1,i2,j1,j2] = get_file_index(Param.IndexRange.first_i,first_j,PairString);
     57    FirstFileName=fullfile_uvmat(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3},...
     58        Param.InputTable{1,5},Param.InputTable{1,4},i1,i2,j1,j2);
     59    return
    4060end
    4161
    42 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%enable waitbar
    43 hseries=guidata(Series.hseries);%handles of the GUI series
    44 WaitbarPos=get(hseries.waitbar_frame,'Position');
    45 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    46 
    47 basename=fullfile(Series.RootPath,Series.RootFile) ;
    48 
    49 %create dir of the new images
    50 [dir_images,namebase]=fileparts(basename);
    51 [path,subdir_ima]=fileparts(dir_images);
    52 curdir=pwd;
    53 cd(path);
    54 mkdir([subdir_ima '_vol']);
    55 cd(curdir);
    56 basename_new=fullfile(path,[subdir_ima '_vol'],namebase);
    57 
    58 % read imadoc
    59 [XmlData,warntext]=imadoc2struct([basename '.xml']);
    60 nbfield1=size(XmlData.Time,1)
    61 nbfield2=size(XmlData.Time,2)
    62 
    63 answer=msgbox_uvmat('INPUT_Y-N','apply image rescaling function levels.m')
    64 test_level=isequal(answer,'Yes')
    65 
    66 %copy the xml file
    67 if exist([basename '.xml'],'file')
    68     copyfile([basename '.xml'],[basename_new '.xml']);% copy the .civ file
    69     t=xmltree([basename_new '.xml']);
    70    
    71     %update information on the first image name in the series
    72     uid_Heading=find(t,'ImaDoc/Heading');
    73     if isempty(uid_Heading)
    74         [t,uid_Heading]=add(t,1,'element','Heading');
    75     end   
    76     uid_ImageName=find(t,'ImaDoc/Heading/ImageName');
    77     ImageName=name_generator(basename_new,num_i1(1),num_j1(1),'.png','_i_j');
    78     [pth,ImageName]=fileparts(ImageName);
    79     ImageName=[ImageName '.png']
    80     if isempty(uid_ImageName)
    81        [t,uid_ImageName]=add(t,uid_Heading,'element','ImageName');
    82     end
    83     uid_value=children(t,uid_ImageName);
    84     if isempty(uid_value)
    85         t=add(t,uid_ImageName,'chardata',ImageName)%indicate  name of the first image, with ;png extension
    86     else
    87         t=set(t,uid_value(1),'value',ImageName)%indicate  name of the first image, with ;png extension
    88     end 
    89     save(t,[basename_new '.xml'])
     62%%%%%%%%%%%% STANDARD PART (DO NOT EDIT) %%%%%%%%%%%%
     63%% read input parameters from an xml file if input is a file name (batch mode)
     64ParamOut=[];
     65RUNHandle=[];
     66WaitbarHandle=[];
     67checkrun=1;
     68if ischar(Param)
     69    Param=xml2struct(Param);% read Param as input file (batch case)
     70    checkrun=0;
     71else% interactive mode in Matlab
     72    hseries=findobj(allchild(0),'Tag','series');
     73    RUNHandle=findobj(hseries,'Tag','RUN');%handle of RUN button in GUI series
     74    WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series
    9075end
    9176
    92 %main loop
    93  vol=[];
    94  for ifile=1:nbfield1*nbfield2
    95      update_waitbar(hseries.waitbar,WaitbarPos,ifile/(nbfield1*nbfield2))
    96      stopstate=get(hseries.RUN,'BusyAction');
    97      if isequal(stopstate,'queue') % enable STOP command
    98          filename=name_generator(basename,ifile-1,1,Series.FileExt,Series.NomType);
    99          num_j=mod(ifile-1,nbfield2)+1;
    100          num_i=floor((ifile-1)/nbfield2)+1;
    101          A=imread(filename);
    102          Atype=class(A);
    103          if test_level
    104              A=levels(A,16);
    105              display(num2str(num_i))
    106          end
    107          vol=[vol;A];%concacene along y
    108          if num_j==nbfield2
    109              filename_new=name_generator(basename_new,num_i,1,'.vol','_i');
    110              imwrite(vol,filename_new,'png','BitDepth',16)% WRITE IN 16 bits: needed for the current version of civ3C3D
    111              display([filename_new 'written (16bits image)'])
    112              vol=[];
    113          end
    114      end
    115  end
     77%% subdirectory for output files
     78SubdirOut=[Param.OutputSubDir Param.OutputDirExt];
     79
     80%% root input file names and nomenclature type (cell arrays with one element)
     81RootPath=Param.InputTable(:,1);
     82RootFile=Param.InputTable(:,3);
     83SubDir=Param.InputTable(:,2);
     84NomType=Param.InputTable(:,4);
     85FileExt=Param.InputTable(:,5);
     86
     87
     88%% get the set of input file names (cell array filecell), and file indices
     89[filecell,i1_series,i2_series,j1_series,j2_series]=get_file_series(Param);
     90% filecell{iview,fileindex}: cell array representing the list of file names
     91%        iview: line in the table corresponding to a given file series
     92%        fileindex: file index within  the file series,
     93% i1_series(iview,ref_j,ref_i)... are the corresponding arrays of indices i1,i2,j1,j2, depending on the input line iview and the two reference indices ref_i,ref_j
     94% i1_series(iview,fileindex) expresses the same indices as a 1D array in file indices
     95nbfield_j=size(i1_series{1},1); %nb of fields for the j index (bursts or volume slices)
     96nbfield_i=size(i1_series{1},2); %nb of fields for the i index
     97nbfield=nbfield_j*nbfield_i; %total number of fields
     98[FileInfo{1},VideoObject{1}]=get_file_info(filecell{1,1});% type of input file
     99FileType{1}=FileInfo{1}.FileType;
     100
     101%% calibration data and timing: read the ImaDoc files
     102[XmlData,NbSlice_calib,time,errormsg]=read_multimadoc(RootPath,SubDir,RootFile,FileExt,i1_series,i2_series,j1_series,j2_series);
     103if ~isempty(errormsg)
     104    disp_uvmat('WARNING',errormsg,checkrun)
     105end
     106
     107%% coordinate transform or other user defined transform
     108transform_fct='';%default fct handle
     109if isfield(Param,'FieldTransform')&&~isempty(Param.FieldTransform.TransformName)
     110        currentdir=pwd;
     111        cd(Param.FieldTransform.TransformPath)
     112        transform_fct=str2func(Param.FieldTransform.TransformName);
     113        cd (currentdir)     
     114end
     115
     116%% main loop
     117for ifile=1:nbfield_i
     118    update_waitbar(WaitbarHandle,ifile/nbfield)
     119    if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue')
     120        disp('program stopped by user')
     121        return
     122    end
     123    for jfile=1:nbfield_j
     124        if ~isempty(j1_series)&&~isequal(j1_series,{[]})
     125            j1=j1_series{1}(jfile,ifile);
     126        end
     127        filename=fullfile_uvmat(RootPath{1},SubDir{1},RootFile{1},FileExt{1},NomType{1},i1_series{1}(jfile,ifile),[],j1)
     128        [Data,tild,errormsg] = read_field(filename,'image');
     129        % transform the input field (e.g; phys) if requested (no transform involving two input fields)
     130        if ~isempty(transform_fct)
     131            Data.ZIndex=jfile;
     132            Data=transform_fct(Data,XmlData{1});
     133        end
     134        if jfile==1
     135            vol=zeros(nbfield_j,size(Data.A,1),size(Data.A,2));
     136            Z=1:nbfield_j;%default Z values
     137        end
     138        vol(jfile,:,:)=double(Data.A);%concacene along y
     139        Z(jfile)=Data.PlaneCoord(3);
     140    end
     141    if ifile==1
     142        npx=size(Data.A,2);
     143        npy=size(Data.A,1);
     144        npz=256;
     145        ind_x=round(npx/2)-10:round(npx/2)+10;%image index at the mid x position
     146        ind_y=round(npy/2)-10:round(npy/2)+10;;%image index at the mid y position
     147       
     148        %write xml calibration file, using the first file
     149            Rangx=Data.Coord_x;
     150            Rangy=Data.Coord_y;
     151            Rangz=[Z(end) Z(1)];
     152   
     153        GeometryCal.CalibrationType='rescale';
     154        GeometryCal.CoordUnit=Data.CoordUnit;
     155        GeometryCal.focal=1;
     156        %scaling along x, y and z
     157        pxcmx=(npx-1)/(Rangx(2)-Rangx(1));
     158        pxcmy=(npy-1)/(Rangy(1)-Rangy(2));
     159        pxcmz=(npz-1)/(Rangz(2)-Rangz(1));
     160        T_x=-pxcmx*Rangx(1)+0.5;
     161        T_y=-pxcmy*Rangy(2)+0.5;
     162        T_z=-pxcmz*Rangz(2)+0.5;
     163        % xml file for x cut
     164        GeometryCal.R=[pxcmx,0,0;0,pxcmz,0;0,0,1];
     165        GeometryCal.Tx_Ty_Tz=[T_x T_z 1];
     166        ImaDoc.GeometryCalib=GeometryCal;
     167        t=struct2xml(ImaDoc);
     168        t=set(t,1,'name','ImaDoc');
     169        save(t,fullfile(RootPath{1},SubdirOut,'cut_x.xml'))
     170                   % xml file for y cut
     171        GeometryCal.R=[pxcmy,0,0;0,pxcmz,0;0,0,1];
     172        GeometryCal.Tx_Ty_Tz=[T_y T_z 1];
     173        ImaDoc.GeometryCalib=GeometryCal;
     174        t=struct2xml(ImaDoc);
     175        t=set(t,1,'name','ImaDoc');
     176        save(t,fullfile(RootPath{1},SubdirOut,'cut_y.xml'))
     177    end
     178    cut_y=squeeze(mean(vol(:,:,ind_x),3));
     179    cut_y=interp1(Z,cut_y,linspace(Z(1),Z(end),npz));
     180    cut_x=squeeze(mean(vol(:,ind_y,:),2));
     181    cut_x=interp1(Z,cut_x,linspace(Z(1),Z(end),npz));
     182   
     183    filename_x=fullfile_uvmat(RootPath{1},SubdirOut,'cut_x','.png','_1',i1_series{1}(jfile,ifile),[],j1);
     184    filename_y=fullfile_uvmat(RootPath{1},SubdirOut,'cut_y','.png','_1',i1_series{1}(jfile,ifile),[],j1);
     185    %  filename_new=name_generator(basename_new,num_i,1,'.vol','_i');
     186    imwrite(uint8(vol),filename_x,'png','BitDepth',8)%
     187    display([filename_x 'written (8bits image)'])
     188    imwrite(uint8(vol),filename_y,'png','BitDepth',8)%
     189    display([filename_y 'written (8bits image)'])
     190end
    116191
    117192
    118193
    119 function C=levels(A,bitdepth)
    120 %whos A;
    121 B=double(A(:,:,1));
    122 windowsize=round(min(size(B,1),size(B,2))/20);
    123 windowsize=floor(windowsize/2)*2+1;
    124 ix=[1/2-windowsize/2:-1/2+windowsize/2];%
    125 %del=np/3;
    126 %fct=exp(-(ix/del).^2);
    127 fct2=cos(ix/(windowsize-1)/2*pi/2);
    128 %Mfiltre=(ones(5,5)/5^2);
    129 %Mfiltre=fct2';
    130 Mfiltre=fct2'*fct2;
    131 Mfiltre=Mfiltre/(sum(sum(Mfiltre)));
    132 
    133 C=filter2(Mfiltre,B);
    134 C(:,1:windowsize)=C(:,windowsize)*ones(1,windowsize);
    135 C(:,end-windowsize+1:end)=C(:,end-windowsize+1)*ones(1,windowsize);
    136 C(1:windowsize,:)=ones(windowsize,1)*C(windowsize,:);
    137 C(end-windowsize+1:end,:)=ones(windowsize,1)*C(end-windowsize,:);
    138 C=tanh(B./(2*C));
    139 [n,c]=hist(reshape(C,1,[]),100);
    140 % figure;plot(c,n);
    141 
    142 [m,i]=max(n);
    143 c_max=c(i);
    144 [dummy,index]=sort(abs(c-c(i)));
    145 n=n(index);
    146 c=c(index);
    147 i_select = find(cumsum(n)<0.95*sum(n));
    148 if isempty(i_select)
    149     i_select = 1:length(c);
    150 end
    151 c_select=c(i_select);
    152 n_select=n(i_select);
    153 cmin=min(c_select);
    154 cmax=max(c_select);
    155 if isequal(bitdepth,16)
    156     C=((C-cmin)/(cmax-cmin))*256*256;
    157     C=uint16(C);
    158 else
    159     C=((C-cmin)/(cmax-cmin))*256;
    160     C=uint8(C);
    161 end
Note: See TracChangeset for help on using the changeset viewer.