source: trunk/src/phys_ima.m @ 1060

Last change on this file since 1060 was 972, checked in by sommeria, 7 years ago

multitif updated

File size: 4.8 KB
Line 
1
2% phys_ima: transform several images in phys coordinates on a common pixel grid
3%------------------------------------------------------------------------
4% OUTPUT:
5% A_out: cell array of oitput images corresponding to the transform of the input images
6% Rangx, Rangy; vectors with two elements defining the phys positions of first and last pixels in each direction
7%  (the same for all the ouput images)
8%
9% INPUT: 
10% A: cell array of input images
11% XmlData: cell array of structures defining the calibration parameters for each image
12% ZIndex: index of the reference plane used to define the phys position in 3D
13
14function [A_out,Rangx,Rangy]=phys_ima(A,XmlData,ZIndex)
15xcorner=[];
16ycorner=[];
17npx=[];
18npy=[];
19dx=ones(1,numel(A));
20dy=ones(1,numel(A));
21if isstruct(XmlData)
22    XmlData={XmlData};
23end
24for icell=1:numel(A)
25    siz=size(A{icell});
26    npx=[npx siz(2)];
27    npy=[npy siz(1)];
28    Calib=XmlData{icell}.GeometryCalib;
29    xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordinates of corners
30    yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5];
31    [xcorner_new,ycorner_new]=phys_XYZ(Calib,xima,yima,ZIndex);%corresponding physical coordinates
32    dx(icell)=(max(xcorner_new)-min(xcorner_new))/(siz(2)-1);
33    dy(icell)=(max(ycorner_new)-min(ycorner_new))/(siz(1)-1);
34    xcorner=[xcorner xcorner_new];
35    ycorner=[ycorner ycorner_new];
36end
37Rangx(1)=min(xcorner);
38Rangx(2)=max(xcorner);
39Rangy(2)=min(ycorner);
40Rangy(1)=max(ycorner);
41test_multi=(max(npx)~=min(npx)) || (max(npy)~=min(npy)); %different image lengths
42% npX=1+round((Rangx(2)-Rangx(1))/min(dx));% nbre of pixels in the new image (use the finest resolution min(dx) in the set of images)
43% npY=1+round((Rangy(1)-Rangy(2))/min(dy));
44
45npX=1+round((Rangx(2)-Rangx(1))/max(dx));% nbre of pixels in the new image (use the largest resolution max(dx) in the set of images)
46npY=1+round((Rangy(1)-Rangy(2))/max(dy));
47
48
49x=linspace(Rangx(1),Rangx(2),npX);
50y=linspace(Rangy(1),Rangy(2),npY);
51[X,Y]=meshgrid(x,y);%grid in physical coordiantes
52%vec_B=[];
53A_out=cell(1,numel(A));
54for icell=1:length(A)
55    Calib=XmlData{icell}.GeometryCalib;
56    % rescaling of the image coordinates without change of the image array
57    if strcmp(Calib.CalibrationType,'rescale') && isequal(Calib,XmlData{1}.GeometryCalib)
58        A_out{icell}=A{icell};%no transform
59        Rangx=[0.5 npx-0.5];%image coordiantes of corners
60        Rangy=[npy-0.5 0.5];
61        [Rangx]=phys_XYZ(Calib,Rangx,[0.5 0.5],ZIndex);%case of translations without rotation and quadratic deformation
62        [xx,Rangy]=phys_XYZ(Calib,[0.5 0.5],Rangy,ZIndex);
63    else
64        % the image needs to be interpolated to the new coordinates
65        zphys=0; %default
66        if isfield(Calib,'SliceCoord') %.Z= index of plane
67            SliceCoord=Calib.SliceCoord(ZIndex,:);
68            zphys=SliceCoord(3); %to generalize for non-parallel planes
69            if isfield(Calib, 'SliceAngle') && ~isequal(Calib.SliceAngle,[0 0 0]) && ~isequal(Calib.SliceAngle(ZIndex,:),[0 0 0])
70                testangle=1;
71                om=norm(Calib.SliceAngle(ZIndex,:));%norm of rotation angle in radians
72                OmAxis=Calib.SliceAngle(ZIndex,:)/om; %unit vector marking the rotation axis
73                cos_om=cos(pi*om/180);
74                sin_om=sin(pi*om/180);
75                coeff=OmAxis(3)*(1-cos_om);
76                norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om;
77                norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om;
78                norm_plane(3)=OmAxis(3)*coeff+cos_om;
79                %Z0=norm_plane*Calib.SliceCoord(ZIndex,:)'/norm_plane(3);
80                Z0=Calib.SliceCoord(ZIndex,3);
81                zphys=Z0-(norm_plane(1)*X-norm_plane(2)*Y)/norm_plane(3);
82            end
83%             if isfield(Calib,'InterfaceCoord') && isfield(Calib,'RefractionIndex')
84%                 H=Calib.InterfaceCoord(3);
85%                 if H>zphys
86%                     zphys=H-(H-zphys)/Calib.RefractionIndex; %corrected z (virtual object)
87%                 end
88%             end
89        end
90        xima=0.5:npx(icell)-0.5;%image coordinates of corners
91        yima=npy(icell)-0.5:-1:0.5;
92        [XIMA_init,YIMA_init]=meshgrid(xima,yima);%grid of initial image in px coordinates
93        [XIMA,YIMA]=px_XYZ(XmlData{icell}.GeometryCalib,X,Y,zphys);% image coordinates for each point in the real
94        testuint8=isa(A{icell},'uint8');
95        testuint16=isa(A{icell},'uint16');
96        if ndims(A{icell})==2 %(B/W images)
97        A_out{icell}=interp2(XIMA_init,YIMA_init,double(A{icell}),XIMA,YIMA);
98         elseif ndims(A{icell})==3     
99             for icolor=1:size(A{icell},3)
100                 A{icell}=double(A{icell});
101                 A_out{icell}(:,:,icolor)=interp2(XIMA_init,YIMA_init,A{icell}(:,:,icolor),XIMA,YIMA);
102             end
103         end
104        if testuint8
105            A_out{icell}=uint8(A_out{icell});
106        end
107        if testuint16
108            A_out{icell}=uint16(A_out{icell});
109        end     
110    end
111end
Note: See TracBrowser for help on using the repository browser.