Ignore:
Timestamp:
Mar 23, 2010, 2:52:49 PM (15 years ago)
Author:
sommeria
Message:

-plot projections on a new specific GUI view_field, to avoid the multiplication of active figures
-introduce a ruler in uvmat (in menu/Tools), to measure distances and angles
-insert automatic detection of points for geometric calibration: tool 'detect_grid'. Four points, deliminating the grid to determine, must be marked with the mouse, as well as the physical grid. Then the points inside are automatically detected.
-reading plane positions in ImaDoc? improved to deal with volume scans

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/geometry_calib.m

    r54 r60  
    4343% Edit the above text to modify the response to help geometry_calib
    4444
    45 % Last Modified by GUIDE v2.5 05-Jan-2010 23:22:04
     45% Last Modified by GUIDE v2.5 23-Mar-2010 06:22:33
    4646
    4747% Begin initialization code - DO NOT edit
    … …  
    7070% PlotHandles: set of handles of the elements contolling the plotting
    7171% parameters on the uvmat interface (obtained by 'get_plot_handle.m')
     72%------------------------------------------------------------------------
    7273function geometry_calib_OpeningFcn(hObject, eventdata, handles, handles_uvmat,pos,inputfile)
    73 
     74%------------------------------------------------------------------------
    7475% Choose default command line output for geometry_calib
    7576handles.output = hObject;
    … …  
    218219
    219220
    220 %----------------------------------------------------
     221%------------------------------------------------------------------------
    221222% executed when closing: set the parent interface button to value 0
    222223function closefcn(gcbo,eventdata,handles_uvmat)
     224%------------------------------------------------------------------------
    223225huvmat=findobj(allchild(0),'Name','uvmat');
    224226if exist('handles_uvmat','var')
    … …  
    229231end
    230232
    231 
    232 % % --- Executes on button press in MenuCoord.
    233 % function MenuCoord_Callback(hObject, eventdata, handles)
    234 
    235 %
    236 % % --- Executes on button press in delete.
    237 % function delete_Callback(hObject, eventdata, handles)
    238 % SetData=get(gcbf,'UserData');%get the interface data
    239 % IndexObj=SetData.IndexObj;
    240 % delete_object(IndexObj);
    241 
    242 
    243 %------------------------------------------------------------------
     233%------------------------------------------------------------------------
    244234% --- Executes on button press in calibrate_lin.
    245235function APPLY_Callback(hObject, eventdata, handles)
    246 %------------------------------------------------------------------
     236%------------------------------------------------------------------------
    247237calib_cell=get(handles.calib_type,'String');
    248238val=get(handles.calib_type,'Value');
    … …  
    280270end
    281271update_imadoc(GeometryCalib,outputfile)
    282 % testappend=0;
    283 % if exist(outputfile,'file');%=1 if the output file already exists, 0 else
    284 %     t=xmltree(outputfile); %read the file
    285 %     backupfile=outputfile;
    286 %     testexist=2;
    287 %     while testexist==2
    288 %         backupfile=[backupfile '~'];% make a backup name by adding  ~ to the xml file name
    289 %         testexist=exist(backupfile,'file');
    290 %     end
    291 %     [success,message]=copyfile(outputfile,backupfile);%make backup   
    292 %     t=xmltree(outputfile); %read the file
    293 %     uid=find(t,'ImaDoc');
    294 %     if ~isequal(uid,1)%if the xml file is not ImaDoc, delete it (after backup)
    295 %         if isequal(success,1)
    296 %             delete(outputfile)
    297 %         else
    298 %             msgbox_uvmat('ERROR',['error in the backup of the existing xml file: ' message])
    299 %             return
    300 %         end
    301 %     else
    302 %         uid_calib=find(t,'ImaDoc/GeometryCalib');
    303 %         testappend=1;
    304 %         if isempty(uid_calib)
    305 %             [t,uid_calib]=add(t,1,'element','GeometryCalib');
    306 %         else %if GeometryCalib already exists, delete its content
    307 %             uid_child=children(t,uid_calib);
    308 %             t=delete(t,uid_child);
    309 % %             testappend=1;
    310 %         end
    311 %     end
    312 % end
    313 % if ~testappend %create a new xml file for calibration data
    314 %     t=xmltree;
    315 %     t=set(t,1,'name','ImaDoc');
    316 %     [t,uid_calib]=add(t,1,'element','GeometryCalib');
    317 % end
    318 % % hgrid=get(handles.REPLICATE,'parent');%read the calibration image source on the interface userdata
    319 % % imagename=get(hgrid,'UserData');
    320 % % if exist(imagename,'file')
    321 % %     GeometryCalib.SourceCalib.ImageCalib=imagename;
    322 % % end
    323 % GeometryCalib.SourceCalib.PointCoord=Object.Coord;
    324 % t=struct2xml(GeometryCalib,t,uid_calib);
    325 % save(t,outputfile);
    326272msgbox_uvmat('CONFIRMATION',{[outputfile ' updated with calibration data'];...
    327273    ['Error rms (along x,y)=' num2str(GeometryCalib.ErrorRms) ' pixels'];...
    … …  
    341287uvmat('RootPath_Callback',hObject,eventdata,hhuvmat); %file input with xml reading  in uvmat
    342288
    343 %------------------------------------------------------------------
     289%------------------------------------------------------------------------
    344290% --- Executes on button press in calibrate_lin.
    345291function REPLICATE_Callback(hObject, eventdata, handles)
    346 %------------------------------------------------------------------
     292%------------------------------------------------------------------------
    347293calib_cell=get(handles.calib_type,'String');
    348294val=get(handles.calib_type,'Value');
    … …  
    388334    testinput=0;
    389335    if isfield(Heading,'SubCampaign') && isequal([filename ext],Heading.SubCampaign)
    390 %         set(hhdataview.RootDirectory,'String',XmlInput)
    391 %         set(hhdataview.SubCampaignTest,'Value',1)
    392336        SubCampaignTest='y';
    393337        testinput=1;
    394338    elseif isfield(Heading,'Campaign') && isequal([filename ext],Heading.Campaign)
    395 %         set(hhdataview.RootDirectory,'String',XmlInput)
    396 %         set(hhdataview.SubCampaignTest,'Value',0)
    397339        testinput=1;
    398340    end
    … …  
    414356    outcome=dataview(XmlInput,SubCampaignTest,GeometryCalib);
    415357end
    416 %     %A COMPLETER
    417 %     dataview('RootDirectory_Callback',hObject,eventdata,hhdataview)
    418 %     ListDevices=get(hhdataview.ListDevices,'String');
    419 %     for ilist=1:length(ListDevices)
    420 %         if isequal(ListDevices{ilist},Device)
    421 %             set(hhdataview.ListDevices,'Value',ilist)
    422 %             dataview('ListDevices_Callback',hObject,eventdata,hhdataview)
    423 %             break
    424 %         end
    425 %     end
    426 
    427 % % hhdataview=guidata(h_dataview);
    428 % CurrentPath=get(hhdataview.RootDirectory,'String');
    429 % ListExperiments=get(hhdataview.ListExperiments,'String');
    430 % Value=get(hhdataview.ListExperiments,'Value');
    431 % if ~isequal(Value,1)
    432 %     ListExperiments=ListExperiments(Value);
    433 % end
    434 % ListDevices=get(hhdataview.ListDevices,'String');
    435 % Value=get(hhdataview.ListDevices,'Value');
    436 % if isequal(Value,1)
    437 %     msgbox_uvmat('ERROR','manually select in the GUI dataview the device being calibrated')
    438 %     return
    439 % else
    440 %     ListDevices=ListDevices(Value);
    441 % end
    442 % ListRecords=get(hhdataview.ListRecords,'String');
    443 % Value=get(hhdataview.ListRecords,'Value');
    444 % if ~isequal(Value,1)
    445 %     ListRecords=ListRecords(Value);
    446 % end
    447 % [ListDevices,ListRecords,ListXml,List]=ListDir(CurrentPath,ListExperiments,ListDevices,ListRecords);
    448 % ListXml=get(hhdataview.ListXml,'String');
    449 % Value=get(hhdataview.ListXml,'Value');
    450 % if isequal(Value,1)
    451 %     msgbox_uvmat('ERROR','you need to select in the GUI dataview the xml files to edit')
    452 %     return
    453 % else
    454 %     ListXml=ListXml(Value);
    455 % end
    456 %
    457 % %update all the selected xml files
    458 % answer=msgbox_uvmat('INPUT_Y-N',[num2str(length(Value)) ' xml files for device ' ListDevices{1} ' will be refreshed with ' calib_type ' calibration data'])
    459 % if ~isequal(answer,'Yes')
    460 %     return
    461 % end
    462 % 'TESTcalib'
    463 % List=DataFiles.List
    464 % for iexp=1:length(List.Experiment)
    465 %     ExpName=List.Experiment{iexp}.name;
    466 %     if isfield(List.Experiment{iexp},'Device')
    467 %         for idevice=1:length(List.Experiment{iexp}.Device)
    468 %             DeviceName=List.Experiment{iexp}.Device{idevice}.name;       
    469 %             if isfield(List.Experiment{iexp}.Device{idevice},'xmlfile')
    470 %                 for ixml=1:length(List.Experiment{iexp}.Device{idevice}.xmlfile)
    471 %                     FileName=List.Experiment{iexp}.Device{idevice}.xmlfile{ixml};
    472 %                     for ilistxml=1:length(ListXml)
    473 %                         if isequal(FileName,ListXml{ilistxml})
    474 %                             set(hhdataview.ListXml,'Value',Value(ilistxml))
    475 %                             drawnow
    476 %                             xmlfullname=fullfile(CurrentPath,ExpName,DeviceName,FileName);
    477 %                             update_imadoc(GeometryCalib,xmlfullname)
    478 %                             break
    479 %                         end
    480 %                     end
    481 %                 end
    482 %              elseif isfield(List.Experiment{iexp}.Device{idevice},'Record')
    483 %                 for irecord=1:length(List.Experiment{iexp}.Device{idevice}.Record)
    484 %                     RecordName=List.Experiment{iexp}.Device{idevice}.Record{irecord}.name;
    485 %                     if isfield(List.Experiment{iexp}.Device{idevice}.Record{irecord},'xmlfile')
    486 %                         for ixml=1:length(List.Experiment{iexp}.Device{idevice}.Record{irecord}.xmlfile)
    487 %                             FileName=List.Experiment{iexp}.Device{idevice}.Record{irecord}.xmlfile{ixml};
    488 %                             for ilistxml=1:length(ListXml)
    489 %                                 if isequal(FileName,ListXml{ilistxml})
    490 %                                     set(hhdataview.ListXml,'Value',Value(ilistxml))
    491 %                                     drawnow
    492 %                                     xmlfullname=fullfile(CurrentPath,ExpName,DeviceName,RecordName,FileName);
    493 %                                     update_imadoc(GeometryCalib,xmlfullname)
    494 %                                     break
    495 %                                 end
    496 %                             end
    497 %                         end
    498 %                     end
    499 %                 end
    500 %             end
    501 %         end
    502 %     end
    503 % end
    504 % set(hhdataview.ListXml,'Value',Value)
    505 
    506 
    507 %-----------------------------------------------------------------
     358
     359%------------------------------------------------------------------------
    508360% determine the parameters for a calibration by an affine function (rescaling and offset, no rotation)
    509361function GeometryCalib=calib_rescale(Coord)
    510 %------------------------------------------------------------------
     362%------------------------------------------------------------------------
    511363 
    512364X=Coord(:,1);
    … …  
    542394GeometryCalib.ErrorMax(2)=max(abs(Ypoints-y_ima));
    543395
    544 
    545 %------------------------------------------------------------------
     396%------------------------------------------------------------------------
    546397% determine the parameters for a calibration by a linear transform matrix (rescale and rotation)
    547398function GeometryCalib=calib_linear(Coord)
    548 %------------------------------------------------------------------
     399%------------------------------------------------------------------------
    549400X=Coord(:,1);
    550401Y=Coord(:,2);
    … …  
    572423GeometryCalib.ErrorMax(2)=max(abs(y1-y_ima));
    573424
    574 
    575 
    576 
    577 %------------------------------------------------------------------
     425%------------------------------------------------------------------------
    578426function GeometryCalib=calib_tsai(Coord)
    579 %------------------------------------------------------------------
     427%------------------------------------------------------------------------
    580428%TSAI
    581429% 'calibration_lin' provides a linear transform on coordinates,
    … …  
    589437    sparam=convert(t);
    590438end
    591 % else
    592 %     %fid = fopen(fullfile(path_UVMAT,'PARAM_WIN.txt'),'r');%open the file with civ binary names
    593 %     xmlfile=fullfile(path_UVMAT,'PARAM_WIN.xml');
    594 %     if exist(xmlfile,'file')
    595 %         t=xmltree(xmlfile);
    596 %         sparam=convert(t);
    597 %     end
    598 % end
    599439if ~isfield(sparam,'GeometryCalib_exe')
    600440    msgbox_uvmat('ERROR',['calibration program <GeometryCalib_exe> undefined in parameter file ' xmlfile])
    … …  
    618458end
    619459calibdat=dlmread('calib.dat');
     460delete('calib.dat')
     461delete('t.txt')
    620462GeometryCalib.CalibrationType='tsai';
    621463GeometryCalib.focal=calibdat(10);
    … …  
    690532
    691533
    692 
     534%------------------------------------------------------------------------
    693535% --- Executes on button press in rotation.
    694536function rotation_Callback(hObject, eventdata, handles)
     537%------------------------------------------------------------------------
    695538angle_rot=(pi/180)*str2num(get(handles.Phi,'String'));
    696539Coord_cell=get(handles.ListCoord,'String');
    … …  
    701544set(handles.YObject,'String',num2str(data.Coord(:,2),4));
    702545
    703 
     546%------------------------------------------------------------------------
    704547function XImage_Callback(hObject, eventdata, handles)
     548%------------------------------------------------------------------------
    705549update_list(hObject, eventdata,handles)
    706550
     551%------------------------------------------------------------------------
    707552function YImage_Callback(hObject, eventdata, handles)
     553%------------------------------------------------------------------------
    708554update_list(hObject, eventdata,handles)
    709555
    … …  
    717563update_list(hObject, eventdata,handles)
    718564
     565%------------------------------------------------------------------------
    719566function update_list(hObject, eventdata, handles)
     567%------------------------------------------------------------------------
    720568str4=get(handles.XImage,'String');
    721569str5=get(handles.YImage,'String');
    … …  
    732580Coord{val}=strline;
    733581set(handles.ListCoord,'String',Coord)
    734 
    735 %--------------------------------------------------------------------
     582%update the plot
     583ListCoord_Callback(hObject, eventdata, handles)
     584
     585%------------------------------------------------------------------------
    736586% --- Executes on selection change in ListCoord.
    737 %--------------------------------------------------------------------
    738587function ListCoord_Callback(hObject, eventdata, handles)
    739 % hObject    handle to ListCoord (see GCBO)
    740 % eventdata  reserved - to be defined in a future version of MATLAB
    741 % handles    structure with handles and user data (see GUIDATA)
    742 
    743 % Hints: contents = get(hObject,'String') returns ListCoord contents as cell array
    744 %        contents{get(hObject,'Value')} returns selected item from ListCoord
    745 %set(handles.edit_append,'Value',2); %set to edit mode
     588%------------------------------------------------------------------------
    746589Coord_cell=get(handles.ListCoord,'String');
    747590val=get(handles.ListCoord,'Value');
    … …  
    788631end
    789632
    790 
    791 %----------------------------------------------------
    792 % --- Executes on button press in rotation_plus.
    793 function rotation_plus_Callback(hObject, eventdata, handles)
    794 Phi=0;
    795 Phi=get(handles.Phi,'String');
    796 if ~isempty(Phi)
    797     Phi=str2num(Phi);
    798 end
    799 rotation(handles,Phi)
    800 
    801 %-------------------------------------------------
    802 % --- Executes on button press in rotation_minus.
    803 function rotation_minus_Callback(hObject, eventdata, handles)
    804 Phi=0;
    805 Phi=get(handles.Phi,'String');
    806 if ~isempty(Phi)
    807     Phi=-str2num(Phi);
    808 end
    809 rotation(handles,Phi)
    810 
    811 
    812 
    813 function O_x_Callback(hObject, eventdata, handles)
    814 
    815 
    816 function O_y_Callback(hObject, eventdata, handles)
    817 
    818 
    819 function O_z_Callback(hObject, eventdata, handles)
    820 
    821 
     633%------------------------------------------------------------------------
    822634% --- Executes on selection change in edit_append.
    823635function edit_append_Callback(hObject, eventdata, handles)
    824 % val=get(handles.PLOT_append,'Value');
    825 % if isequal(val,2); %append mode
    826 %     %appeler mouse
    827 % end
     636%------------------------------------------------------------------------
    828637choice=get(handles.edit_append,'Value');
    829638if choice==1
    … …  
    839648
    840649
    841 %A REVOIR
    842 % if choice==2
    843 %     %display image with px coordinates
    844 %     hrootpath=findobj(huvmat,'Tag','RootPath');
    845 %     hrootfile=findobj(huvmat,'Tag','RootFile');
    846 %     RootPath='';
    847 %     RootFile='';
    848 % %     if ~isempty(hrootpath)& ~isempty(hrootfile)
    849 %         testhandle=1;
    850 %         RootPath=get(hrootpath,'String');
    851 %         RootFile=get(hrootfile,'String');
    852 % %         filebase=fullfile(RootPath,RootFile);
    853 % %         outputfile=[filebase '.xml'];
    854 %         Indices=get(findobj(huvmat,'Tag','FileIndex'),'String');
    855 %         Ext=get(findobj(huvmat,'Tag','FileExt'),'String');
    856 %         imagename=[fullfile(RootPath,RootFile) Indices Ext];
    857 %         % input.menu_coord=1;
    858 %          h_menu_coord=findobj(huvmat,'Tag','menu_coord');
    859 %         set(h_menu_coord,'Value',3)
    860 %         huvmat=uvmat(imagename);%open uvmat, set phys coord (Value 1)
    861 %     
    862 % %     end
    863 % end
    864650   
    865651function NEW_Callback(hObject, eventdata, handles)
    … …  
    900686 
    901687
    902 
    903 %'key_press_fcn:' function activated when a key is pressed on the keyboard
    904 %-----------------------------------
     688%------------------------------------------------------------------------
     689% --- 'key_press_fcn:' function activated when a key is pressed on the keyboard
    905690function key_press_fcn(hObject,eventdata,handles)
     691%------------------------------------------------------------------------
    906692hh=get(hObject,'parent');
    907693xx=double(get(hh,'CurrentCharacter')); %get the keyboard character
    … …  
    930716end
    931717
    932 
     718%------------------------------------------------------------------------
    933719% --- Executes on button press in append_point.
    934720function append_point_Callback(hObject, eventdata, handles)
    935 
     721%------------------------------------------------------------------------
    936722       Coord=get(handles.ListCoord,'String');
    937723       val=length(Coord);
    … …  
    943729       set(handles.ListCoord,'Value',val+1)
    944730
    945 
    946 % --------------------------------------------------------------------
     731%------------------------------------------------------------------------
    947732function MenuOpen_Callback(hObject, eventdata, handles)
     733%------------------------------------------------------------------------
    948734%get the object file
    949735huvmat=findobj(allchild(0),'Name','uvmat');
    … …  
    972758
    973759
    974 % --------------------------------------------------------------------
    975 function Untitled_3_Callback(hObject, eventdata, handles)
    976 % hObject    handle to Untitled_3 (see GCBO)
    977 % eventdata  reserved - to be defined in a future version of MATLAB
    978 % handles    structure with handles and user data (see GUIDATA)
    979 
    980 
    981 % --------------------------------------------------------------------
     760%------------------------------------------------------------------------
    982761function MenuPlot_Callback(hObject, eventdata, handles)
    983 
     762%------------------------------------------------------------------------
    984763huvmat=findobj(allchild(0),'Name','uvmat');%find the current uvmat interface handle
    985764UvData=get(huvmat,'UserData');%Data associated to the current uvmat interface
    986765hhuvmat=guidata(huvmat); %handles of GUI elements in uvmat
    987766hplot=findobj(huvmat,'Tag','axes3');%main plotting axis of uvmat
    988 h_menu_coord=findobj(huvmat,'Tag','menu_coord');
     767h_menu_coord=findobj(huvmat,'Tag','transform_fct');
    989768menu=get(h_menu_coord,'String');
    990769choice=get(h_menu_coord,'Value');
    … …  
    1033812    Tinput=CalibData.grid;
    1034813end
    1035 T=create_grid(Tinput);%display the GUI create_grid
    1036 CalibData.grid=T;
     814[T,CalibData.grid]=create_grid(grid_input);%display the GUI create_grid
    1037815set(handles.figure1,'UserData',CalibData)
    1038816
    … …  
    1134912
    1135913
     914% --------------------------------------------------------------------
     915function MenuDetectGrid_Callback(hObject, eventdata, handles)
     916
     917CalibData=get(handles.figure1,'UserData');
     918grid_input=[];%default
     919if isfield(CalibData,'grid')
     920    grid_input=CalibData.grid;%retrieve the previously used grid
     921end
     922[T,CalibData.grid]=create_grid(grid_input);%display the GUI create_grid
     923set(handles.figure1,'UserData',CalibData)%store the phys grid for later use
     924
     925%read the four last point coordiantes in pixels
     926Coord_cell=get(handles.ListCoord,'String');%read list of coordiantes on geometry_calib
     927data=read_geometry_calib(Coord_cell);
     928nbpoints=size(data.Coord,1); %nbre of calibration points
     929if nbpoints<4
     930    msgbox_uvmat('ERROR','four points must be selected by the mouse to delimitate the detection area')
     931end
     932corners_X=(data.Coord(end-3:end,4)); %pixel absissa of the four corners
     933corners_Y=(data.Coord(end-3:end,5));
     934
     935%read the current image
     936huvmat=findobj(allchild(0),'Name','uvmat');
     937UvData=get(huvmat,'UserData');
     938A=UvData.Field.A;
     939npxy=size(A);
     940%linear transform on the current image
     941X=[CalibData.grid.x_0 CalibData.grid.x_1 CalibData.grid.x_0 CalibData.grid.x_1]';%corner absissa in the rectified image
     942Y=[CalibData.grid.y_0 CalibData.grid.y_0 CalibData.grid.y_1 CalibData.grid.y_1]';%corner absissa in the rectified image
     943XY_mat=[ones(size(X)) X Y];
     944a_X1=XY_mat\corners_X; %transformation matrix for X
     945x1=XY_mat*a_X1;%reconstruction
     946err_X1=max(abs(x1-corners_X))%error
     947a_Y1=XY_mat\corners_Y;%transformation matrix for X
     948y1=XY_mat*a_Y1;
     949err_Y1=max(abs(y1-corners_Y))%error
     950GeometryCalib.CalibrationType='linear';
     951GeometryCalib.CoordUnit=[];% default value, to be updated by the calling function
     952GeometryCalib.f=1;
     953GeometryCalib.dpx=1;
     954GeometryCalib.dpy=1;
     955GeometryCalib.sx=1;
     956GeometryCalib.Cx=0;
     957GeometryCalib.Cy=0;
     958GeometryCalib.kappa1=0;
     959GeometryCalib.Tx=a_X1(1);
     960GeometryCalib.Ty=a_Y1(1);
     961GeometryCalib.Tz=1;
     962GeometryCalib.R=[a_X1(2),a_X1(3),0;a_Y1(2),a_Y1(3),0;0,0,1];
     963[Amod,Rangx,Rangy]=phys_Ima(A-min(min(A)),GeometryCalib,0);
     964Amod=double(Amod);
     965%figure(12)
     966%Amax=max(max(Amod))
     967%image(Rangx,Rangy,uint8(255*Amod/Amax))
     968ind_range=10;% range of search of image ma around each point obtained by linear interpolation from the marked points
     969nbpoints=size(T,1);
     970for ipoint=1:nbpoints
     971    Dx=(Rangx(2)-Rangx(1))/(npxy(2)-1); %x mesh in real space
     972    Dy=(Rangy(2)-Rangy(1))/(npxy(1)-1); %y mesh in real space
     973    i0=1+round((T(ipoint,1)-Rangx(1))/Dx);%round(Xpx(ipoint));
     974    j0=1+round((T(ipoint,2)-Rangy(1))/Dy);%round(Xpx(ipoint));
     975    Asub=Amod(j0-ind_range:j0+ind_range,i0-ind_range:i0+ind_range);
     976    x_profile=sum(Asub,1);
     977    y_profile=sum(Asub,2);
     978    [Amax,ind_x_max]=max(x_profile);
     979    [Amax,ind_y_max]=max(y_profile);
     980    Delta(ipoint,1)=(ind_x_max-ind_range-1)*Dx;%shift from the initial guess
     981    Delta(ipoint,2)=(ind_y_max-ind_range-1)*Dy;
     982end
     983Tmod=T(:,(1:2))+Delta;
     984[Xpx,Ypx]=px_XYZ(GeometryCalib,Tmod(:,1),Tmod(:,2));
     985for i=1:nbpoints
     986     Coord{i,1}=num2str(T(i,1),4);%display coordiantes with 4 digits
     987     Coord{i,2}=num2str(T(i,2),4);%display coordiantes with 4 digits
     988     Coord{i,3}='0';
     989     Coord{i,4}=num2str(Xpx(i),4);%display coordiantes with 4 digi
     990     Coord{i,5}=num2str(Ypx(i),4);%display coordiantes with 4 digi
     991end
     992Tabchar=cell2tab(Coord,'    |    ');
     993set(handles.ListCoord,'Value',1)
     994set(handles.ListCoord,'String',Tabchar)
     995
     996
     997%%%%%%%%%%%%%%%%%%%%
     998function [A_out,Rangx,Rangy]=phys_Ima(A,Calib,ZIndex)
     999xcorner=[];
     1000ycorner=[];
     1001npx=[];
     1002npy=[];
     1003siz=size(A);
     1004npx=[npx siz(2)];
     1005npy=[npy siz(1)];
     1006xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordiantes of corners
     1007yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
     1008[xcorner,ycorner]=phys_XYZ(Calib,xima,yima,ZIndex);%corresponding physical coordinates
     1009Rangx(1)=min(xcorner);
     1010Rangx(2)=max(xcorner);
     1011Rangy(2)=min(ycorner);
     1012Rangy(1)=max(ycorner);
     1013test_multi=(max(npx)~=min(npx)) | (max(npy)~=min(npy));
     1014npx=max(npx);
     1015npy=max(npy);
     1016x=linspace(Rangx(1),Rangx(2),npx);
     1017y=linspace(Rangy(1),Rangy(2),npy);
     1018[X,Y]=meshgrid(x,y);%grid in physical coordiantes
     1019vec_B=[];
     1020
     1021zphys=0; %default
     1022if isfield(Calib,'SliceCoord') %.Z= index of plane
     1023   SliceCoord=Calib.SliceCoord(ZIndex,:);
     1024   zphys=SliceCoord(3); %to generalize for non-parallel planes
     1025end
     1026[XIMA,YIMA]=px_XYZ(Calib,X,Y,zphys);%corresponding image indices for each point in the real space grid
     1027XIMA=reshape(round(XIMA),1,npx*npy);%indices reorganized in 'line'
     1028YIMA=reshape(round(YIMA),1,npx*npy);
     1029flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
     1030testuint8=isa(A,'uint8');
     1031testuint16=isa(A,'uint16');
     1032if numel(siz)==2 %(B/W images)
     1033    vec_A=reshape(A,1,npx*npy);%put the original image in line
     1034    ind_in=find(flagin);
     1035    ind_out=find(~flagin);
     1036    ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
     1037    ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
     1038    vec_B(ind_in)=vec_A(ICOMB);
     1039    vec_B(ind_out)=zeros(size(ind_out));
     1040    A_out=reshape(vec_B,npy,npx);%new image in real coordinates
     1041elseif numel(siz)==3     
     1042    for icolor=1:siz(3)
     1043        vec_A=reshape(A{icell}(:,:,icolor),1,npx*npy);%put the original image in line
     1044        ind_in=find(flagin);
     1045        ind_out=find(~flagin);
     1046        ICOMB=((XIMA-1)*npy+(npy+1-YIMA));
     1047        ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
     1048        vec_B(ind_in)=vec_A(ICOMB);
     1049        vec_B(ind_out)=zeros(size(ind_out));
     1050        A_out(:,:,icolor)=reshape(vec_B,npy,npx);%new image in real coordinates
     1051    end
     1052end
     1053if testuint8
     1054    A_out=uint8(A_out);
     1055end
     1056if testuint16
     1057    A_out=uint16(A_out);
     1058end
     1059
     1060%INPUT:
     1061%Z: index of plane
     1062function [Xphys,Yphys,Zphys]=phys_XYZ(Calib,X,Y,Z)
     1063if exist('Z','var')& isequal(Z,round(Z))& Z>0 & isfield(Calib,'SliceCoord')&length(Calib.SliceCoord)>=Z
     1064    Zindex=Z;
     1065    Zphys=Calib.SliceCoord(Zindex,3);%GENERALISER AUX CAS AVEC ANGLE
     1066else
     1067%     if exist('Z','var')
     1068%         Zphys=Z;
     1069%     else
     1070        Zphys=0;
     1071%     end
     1072end
     1073if ~exist('X','var')||~exist('Y','var')
     1074    Xphys=[];
     1075    Yphys=[];%default
     1076    return
     1077end
     1078Xphys=X;%default
     1079Yphys=Y;
     1080%image transform
     1081if isfield(Calib,'R')
     1082    R=(Calib.R)';
     1083    Dx=R(5)*R(7)-R(4)*R(8);
     1084    Dy=R(1)*R(8)-R(2)*R(7);
     1085    D0=Calib.f*(R(2)*R(4)-R(1)*R(5));
     1086    Z11=R(6)*R(8)-R(5)*R(9);
     1087    Z12=R(2)*R(9)-R(3)*R(8); 
     1088    Z21=R(4)*R(9)-R(6)*R(7);
     1089    Z22=R(3)*R(7)-R(1)*R(9);
     1090    Zx0=R(3)*R(5)-R(2)*R(6);
     1091    Zy0=R(1)*R(6)-R(3)*R(4);
     1092    A11=R(8)*Calib.Ty-R(5)*Calib.Tz+Z11*Zphys;
     1093    A12=R(2)*Calib.Tz-R(8)*Calib.Tx+Z12*Zphys;
     1094    A21=-R(7)*Calib.Ty+R(4)*Calib.Tz+Z21*Zphys;
     1095    A22=-R(1)*Calib.Tz+R(7)*Calib.Tx+Z11*Zphys;
     1096    X0=Calib.f*(R(5)*Calib.Tx-R(2)*Calib.Ty+Zx0*Zphys);
     1097    Y0=Calib.f*(-R(4)*Calib.Tx+R(1)*Calib.Ty+Zy0*Zphys);
     1098        %px to camera:
     1099    Xd=(Calib.dpx/Calib.sx)*(X-Calib.Cx); % sensor coordinates
     1100    Yd=Calib.dpy*(Y-Calib.Cy);
     1101    dist_fact=1+Calib.kappa1*(Xd.*Xd+Yd.*Yd); %distortion factor
     1102    Xu=dist_fact.*Xd;%undistorted sensor coordinates
     1103    Yu=dist_fact.*Yd;
     1104    denom=Dx*Xu+Dy*Yu+D0;
     1105    % denom2=denom.*denom;
     1106    Xphys=(A11.*Xu+A12.*Yu+X0)./denom;%world coordinates
     1107    Yphys=(A21.*Xu+A22.*Yu+Y0)./denom;
     1108end
Note: See TracChangeset for help on using the changeset viewer.