source: trunk/src/phys_ima.m @ 1134

Last change on this file since 1134 was 1129, checked in by sommeria, 7 months ago

LIF calibration update and extrct_multitif

File size: 4.6 KB
RevLine 
[892]1
[849]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
[1129]13% resolution_factor: factor to control the number of pixels of the new image, default value =1 : same nbre of pixels as the original
14function [A_out,Rangx,Rangy]=phys_ima(A,XmlData,ZIndex,resolution_factor)
[849]15xcorner=[];
16ycorner=[];
17npx=[];
18npy=[];
19dx=ones(1,numel(A));
20dy=ones(1,numel(A));
21if isstruct(XmlData)
22    XmlData={XmlData};
23end
[1129]24if ~exist('resolution_factor','var')
25resolution_factor=1;
26end
[1078]27
[849]28for icell=1:numel(A)
29    siz=size(A{icell});
30    npx=[npx siz(2)];
31    npy=[npy siz(1)];
[1112]32    Calib{icell}=XmlData{icell}.GeometryCalib;
33    Slice{icell}=Calib{icell};
34    if isfield(XmlData{icell},'Slice')
35    Slice{icell}=XmlData{icell}.Slice;
36    end
[1078]37    coord_x=[0.5 siz(2)-0.5];
38    coord_y=[0.5 siz(1)-0.5];
39    x_edge=[linspace(coord_x(1),coord_x(end),npx(icell)) coord_x(end)*ones(1,npy(icell))...
40        linspace(coord_x(end),coord_x(1),npx(icell)) coord_x(1)*ones(1,npy(icell))];%x coordinates of the image edge(four sides)
41    y_edge=[coord_y(1)*ones(1,npx(icell)) linspace(coord_y(1),coord_y(end),npy(icell))...
42        coord_y(end)*ones(1,npx(icell)) linspace(coord_y(end),coord_y(1),npy(icell))];%y coordinates of the image edge(four sides)
[1112]43    [xcorner_new,ycorner_new]=phys_XYZ(Calib{icell},Slice{icell},x_edge,y_edge,ZIndex);%corresponding physical coordinates
[849]44    dx(icell)=(max(xcorner_new)-min(xcorner_new))/(siz(2)-1);
45    dy(icell)=(max(ycorner_new)-min(ycorner_new))/(siz(1)-1);
46    xcorner=[xcorner xcorner_new];
47    ycorner=[ycorner ycorner_new];
48end
49Rangx(1)=min(xcorner);
50Rangx(2)=max(xcorner);
51Rangy(2)=min(ycorner);
52Rangy(1)=max(ycorner);
53test_multi=(max(npx)~=min(npx)) || (max(npy)~=min(npy)); %different image lengths
[922]54
[1129]55npX=1+round(resolution_factor*(Rangx(2)-Rangx(1))/max(dx));% nbre of pixels in the new image (use the largest resolution max(dx) in the set of images)
56npY=1+round(resolution_factor*(Rangy(1)-Rangy(2))/max(dy));
[922]57
[849]58x=linspace(Rangx(1),Rangx(2),npX);
59y=linspace(Rangy(1),Rangy(2),npY);
[1078]60[X,Y]=meshgrid(x,y);%grid in physical coordinates
[849]61A_out=cell(1,numel(A));
[1078]62
[1112]63for icell=1:numel(A)
[849]64    % rescaling of the image coordinates without change of the image array
[1112]65    if strcmp(Calib{icell}.CalibrationType,'rescale') && isequal(Calib,XmlData{1}.GeometryCalib)
[849]66        A_out{icell}=A{icell};%no transform
67        Rangx=[0.5 npx-0.5];%image coordiantes of corners
68        Rangy=[npy-0.5 0.5];
[1112]69        [Rangx]=phys_XYZ(Calib{icell},[],Rangx,[0.5 0.5],ZIndex);%case of translations without rotation and quadratic deformation
70        [~,Rangy]=phys_XYZ(Calib{icell},[],[0.5 0.5],Rangy,ZIndex);
[972]71    else
[849]72        % the image needs to be interpolated to the new coordinates
[1078]73        Z=0; %default
[1112]74        if isfield(Slice{icell},'SliceCoord')&& size(Slice{icell}.SliceCoord,1)>=ZIndex %.Z= index of plane
75            SliceCoord=Slice{icell}.SliceCoord(ZIndex,:);
[1078]76            Z=SliceCoord(3);
[1112]77            if isfield(Slice{icell}, 'SliceAngle') && size(Slice{icell}.SliceAngle,1)>=ZIndex && ~isequal(Slice{icell}.SliceAngle(ZIndex,:),[0 0 0])
78                norm_plane=angle2normal(Slice{icell}.SliceAngle(ZIndex,:));
[1078]79                Z=Z-(norm_plane(1)*(X-SliceCoord(1))+norm_plane(2)*(Y-SliceCoord(2)))/norm_plane(3);
[972]80            end
[849]81        end
[922]82        xima=0.5:npx(icell)-0.5;%image coordinates of corners
83        yima=npy(icell)-0.5:-1:0.5;
[849]84        [XIMA_init,YIMA_init]=meshgrid(xima,yima);%grid of initial image in px coordinates
[1112]85        [XIMA,YIMA]=px_XYZ(Calib{icell},Slice{icell},X,Y,Z);% image coordinates for each point in the real
[849]86        testuint8=isa(A{icell},'uint8');
87        testuint16=isa(A{icell},'uint16');
[1078]88        if ismatrix(A{icell}) %(B/W images)
89            A_out{icell}=interp2(XIMA_init,YIMA_init,double(A{icell}),XIMA,YIMA);
90        elseif ndims(A{icell})==3
91            for icolor=1:size(A{icell},3)
92                A{icell}=double(A{icell});
93                A_out{icell}(:,:,icolor)=interp2(XIMA_init,YIMA_init,A{icell}(:,:,icolor),XIMA,YIMA);
94            end
95        end
[849]96        if testuint8
97            A_out{icell}=uint8(A_out{icell});
98        end
99        if testuint16
100            A_out{icell}=uint16(A_out{icell});
[1078]101        end
[849]102    end
103end
Note: See TracBrowser for help on using the repository browser.