Index: trunk/src/geometry_calib.m
===================================================================
--- trunk/src/geometry_calib.m	(revision 60)
+++ trunk/src/geometry_calib.m	(revision 61)
@@ -43,5 +43,5 @@
 % Edit the above text to modify the response to help geometry_calib
 
-% Last Modified by GUIDE v2.5 23-Mar-2010 06:22:33
+% Last Modified by GUIDE v2.5 23-Mar-2010 16:19:01
 
 % Begin initialization code - DO NOT edit
@@ -148,5 +148,5 @@
 %case of calibration (ImaDoc) input file
 % hcalib=get(handles.calib_type,'parent');
-CalibData=get(handles.figure1,'UserData');
+CalibData=get(handles.geometry_calib,'UserData');
 CalibData.XmlInput=fileinput;
 if isfield(s,'Heading')
@@ -154,5 +154,5 @@
 end
 
-set(handles.figure1,'UserData',CalibData);%store the heading in the interface 'UserData'
+set(handles.geometry_calib,'UserData',CalibData);%store the heading in the interface 'UserData'
 if isfield(s,'GeometryCalib')
     Calib=s.GeometryCalib;
@@ -312,5 +312,5 @@
     delete(h_dataview)
 end
-CalibData=get(handles.figure1,'UserData');%read the calibration image source on the interface userdata
+CalibData=get(handles.geometry_calib,'UserData');%read the calibration image source on the interface userdata
 
 if isfield(CalibData,'XmlInput')
@@ -791,4 +791,6 @@
     set(hh,'YData',ObjectData.Coord(:,2))
 end
+pause(.1)
+figure(handles.geometry_calib)
 
 % --------------------------------------------------------------------
@@ -807,5 +809,5 @@
 %------------------------------------------------------------------------
 %hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib
-CalibData=get(handles.figure1,'UserData');
+CalibData=get(handles.geometry_calib,'UserData');
 Tinput=[];%default
 if isfield(CalibData,'grid')
@@ -813,5 +815,5 @@
 end
 [T,CalibData.grid]=create_grid(grid_input);%display the GUI create_grid
-set(handles.figure1,'UserData',CalibData)
+set(handles.geometry_calib,'UserData',CalibData)
 
 %grid in phys space
@@ -844,5 +846,5 @@
 %-----------------------------------------------------------------------
 %hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib
-CalibData=get(handles.figure1,'UserData');
+CalibData=get(handles.geometry_calib,'UserData');
 Tinput=[];%default
 if isfield(CalibData,'translate')
@@ -851,5 +853,5 @@
 T=translate_points(Tinput);%display translate_points GUI and get shift parameters 
 CalibData.translate=T;
-set(handles.figure1,'UserData',CalibData)
+set(handles.geometry_calib,'UserData',CalibData)
 %translation
 Coord_cell=get(handles.ListCoord,'String');
@@ -872,5 +874,5 @@
 function MenuRotatePoints_Callback(hObject, eventdata, handles)
 %hcalib=get(handles.calib_type,'parent');%handles of the GUI geometry_calib
-CalibData=get(handles.figure1,'UserData');
+CalibData=get(handles.geometry_calib,'UserData');
 Tinput=[];%default
 if isfield(CalibData,'rotate')
@@ -879,5 +881,5 @@
 T=rotate_points(Tinput);%display translate_points GUI and get shift parameters 
 CalibData.rotate=T;
-set(handles.figure1,'UserData',CalibData)
+set(handles.geometry_calib,'UserData',CalibData)
 %-----------------------------------------------------
 %rotation
@@ -915,5 +917,5 @@
 function MenuDetectGrid_Callback(hObject, eventdata, handles)
 
-CalibData=get(handles.figure1,'UserData');
+CalibData=get(handles.geometry_calib,'UserData');
 grid_input=[];%default
 if isfield(CalibData,'grid')
@@ -921,5 +923,5 @@
 end
 [T,CalibData.grid]=create_grid(grid_input);%display the GUI create_grid 
-set(handles.figure1,'UserData',CalibData)%store the phys grid for later use
+set(handles.geometry_calib,'UserData',CalibData)%store the phys grid for later use
 
 %read the four last point coordiantes in pixels
@@ -978,6 +980,17 @@
     [Amax,ind_x_max]=max(x_profile);
     [Amax,ind_y_max]=max(y_profile);
-    Delta(ipoint,1)=(ind_x_max-ind_range-1)*Dx;%shift from the initial guess
-    Delta(ipoint,2)=(ind_y_max-ind_range-1)*Dy;
+    %sub-pixel improvement using moments
+    x_shift=0;
+    y_shift=0;
+    if ind_x_max+2<=2*ind_range+1 && ind_x_max-2>=1
+        Atop=x_profile(ind_x_max-2:ind_x_max+2);
+        x_shift=sum(Atop.*[-2 -1 0 1 2])/sum(Atop);
+    end
+    if ind_y_max+2<=2*ind_range+1 && ind_y_max-2>=1
+        Atop=y_profile(ind_y_max-2:ind_y_max+2);
+        y_shift=sum(Atop.*[-2 -1 0 1 2]')/sum(Atop);
+    end
+    Delta(ipoint,1)=(x_shift+ind_x_max-ind_range-1)*Dx;%shift from the initial guess
+    Delta(ipoint,2)=(y_shift+ind_y_max-ind_range-1)*Dy;
 end
 Tmod=T(:,(1:2))+Delta;
