- Timestamp:
- May 30, 2018, 7:30:21 PM (6 years ago)
- Location:
- trunk/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/transform_field/phys_polar.m
r1027 r1046 95 95 %parameters for polar coordinates (taken from the calibration data of the first field) 96 96 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 97 XmlData.PolarReferenceRadius= 450;98 XmlData.PolarReferenceAngle= 450*pi/2;97 XmlData.PolarReferenceRadius=0;%450; 98 XmlData.PolarReferenceAngle=0;%450*pi/2; 99 99 origin_xy=[0 0];%center for the polar coordinates in the original x,y coordinates 100 100 radius_offset=0;%reference radius used to offset the radial coordinate r 101 101 angle_offset=0; %reference angle used as new origin of the polar angle (= axis Ox by default) 102 angle_scale=180/pi; 102 103 if isfield(XmlData,'TransformInput') 103 104 if isfield(XmlData.TransformInput,'PolarCentre') && isnumeric(XmlData.TransformInput.PolarCentre) … … 249 250 npx=[]; 250 251 npy=[]; 252 NbPoints=20; % nbre of points used to determine the image edge 251 253 for icell=1:length(A) 252 254 siz=size(A{icell}); … … 255 257 zphys=0; %default 256 258 if isfield(CalibIn{icell},'SliceCoord') %.Z= index of plane 259 if ZIndex==0 260 ZIndex=1; 261 end 257 262 SliceCoord=CalibIn{icell}.SliceCoord(ZIndex,:); 258 263 zphys=SliceCoord(3); %to generalize for non-parallel planes 259 264 end 260 xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordiantes of corners 261 yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5]; 265 % xima=[0.5 siz(2)-0.5 0.5 siz(2)-0.5];%image coordinates of corners 266 % yima=[0.5 0.5 siz(1)-0.5 siz(1)-0.5]; 267 edge_x=linspace(0.5,siz(1)-0.5,NbPoints); 268 edge_y=linspace(0.5,siz(2)-0.5,NbPoints); 269 xima=[edge_y (siz(2)-0.5)*ones(1,NbPoints) edge_y 0.5*ones(1,NbPoints)];%image coordinates of corners 270 yima=[0.5*ones(1,NbPoints) edge_x (siz(1)-0.5)*ones(1,NbPoints) edge_x];%image coordinates of corners 262 271 [xcorner_new,ycorner_new]=phys_XYZ(CalibIn{icell},xima,yima,ZIndex);%corresponding physical coordinates 263 272 %transform the corner coordinates into polar ones -
trunk/src/uvmat.m
r1045 r1046 1370 1370 end 1371 1371 1372 %% read the lines previously used for LIF calibration if available, and display them 1373 if isfield(XmlData,'LIFCalib') 1374 if isfield(XmlData.LIFCalib,'Ray1Coord') && isfield(XmlData.LIFCalib,'Ray2Coord') && isfield(XmlData.LIFCalib,'RefLineCoord') 1375 data.Name='Ray1'; 1376 data.Type='line'; 1377 data.ProjMode='none';%default 1378 data.Coord=XmlData.LIFCalib.Ray1Coord; 1379 UvData.ProjObject{2}=data; 1380 UvData.ProjObject{2}.DisplayHandle.uvmat=plot_object(UvData.ProjObject{2},UvData.ProjObject{1},handles.PlotAxes,'b'); 1381 data.Name='Ray2'; 1382 data.Coord=XmlData.LIFCalib.Ray2Coord; 1383 UvData.ProjObject{3}=data; 1384 UvData.ProjObject{3}.DisplayHandle.uvmat=plot_object(UvData.ProjObject{3},UvData.ProjObject{1},handles.PlotAxes,'b'); 1385 data.Name='RefLine'; 1386 data.Coord=XmlData.LIFCalib.RefLineCoord; 1387 UvData.ProjObject{4}=data; 1388 UvData.ProjObject{4}.DisplayHandle.uvmat=plot_object(UvData.ProjObject{4},UvData.ProjObject{1},handles.PlotAxes,'b'); 1389 end 1390 if isfield(XmlData.LIFCalib,'MaskPolygon') 1391 data.Name='MaskPolygon'; 1392 data.Coord=XmlData.LIFCalib.MaskPolygonCoord; 1393 UvData.ProjObject{5}=data; 1394 UvData.ProjObject{5}.DisplayHandle.uvmat=plot_object(UvData.ProjObject{5},UvData.ProjObject{1},handles.PlotAxes,'b'); 1395 end 1396 end 1397 set(handles.ListObject,'String',{'plane';'Ray1';'Ray2';'RefLine';'MaskPolygon'}) 1398 set(handles.ListObject,'Value',1) 1399 set(handles.uvmat,'UserData',UvData); 1400 1372 1401 %% read lines currently drawn 1373 ListObj=UvData.ProjObject; 1402 ListObj=UvData.ProjObject;% get the current list of projection objects 1374 1403 select=zeros(1,numel(ListObj)); 1404 select_mask=zeros(1,numel(ListObj)); 1375 1405 for iobj=1:numel(ListObj); 1376 1406 if isfield(ListObj{iobj},'Type') && strcmp(ListObj{iobj}.Type,'line') 1377 select(iobj)=1; 1407 select(iobj)=1;% select the lines among the projection objects 1408 end 1409 if isfield(ListObj{iobj},'Type') && strcmp(ListObj{iobj}.Type,'polygon') 1410 select_mask(iobj)=1;% select the mask among the projection objects 1378 1411 end 1379 1412 end 1380 1413 val=find(select); 1381 1414 if numel(val)<2 1382 msgbox_uvmat('ERROR','light rays must be defined by at least two lines created by Projection object/line in the menu bar'); 1415 msgbox_uvmat('ERROR',{'light rays must be defined by at least two lines created by Projection object/line in the menu bar'; ... 1416 'use a third line to get a reference luminosity profile accross the illumination beam'}); 1383 1417 return 1384 1418 else 1385 set(handles.ListObject,'Value',val);% show the selected lines on the list1386 1419 ObjectData=UvData.ProjObject(val); 1387 1420 for iobj=1:length(ObjectData) … … 1391 1424 yB(iobj)=ObjectData{iobj}.Coord(2,2); 1392 1425 end 1426 end 1427 index_mask=find(select_mask); 1428 if numel(index_mask)>=1 1429 MaskData=UvData.ProjObject(index_mask); 1430 for iobj=1:length(MaskData) 1431 1432 end 1433 end 1434 1435 %% set the image offset 1436 blackoffset=0; 1437 if isfield(XmlData,'LIFCalib')&& isfield(XmlData.LIFCalib,'BlackOffset') 1438 blackoffset=XmlData.LIFCalib.BlackOffset ;% image value for black background, to be determined by taking images with a cover on the objective lens 1439 end 1440 answer=msgbox_uvmat('INPUT_TXT','camera offset value in the absence of illumination:', num2str(blackoffset)); 1441 if strcmp(answer,'Cancel') 1442 return 1443 else 1444 XmlData.LIFCalib.BlackOffset=answer ;% image value for black background, to be determined by taking images with a cover on the objective lens 1393 1445 end 1394 1446 … … 1405 1457 x0=((x3-x4)*(x1*y2-y1*x2)-(x1-x2)*(x3*y4-y3*x4))/D; 1406 1458 y0=((y3-y4)*(x1*y2-y1*x2)-(y1-y2)*(x3*y4-y3*x4))/D; 1407 XmlData.Illumination.Origin=[x0 y0]; 1408 XmlData.PolarCentre=[x0 y0]; 1459 XmlData.LIFCalib.LightOrigin=[x0 y0];% origin of the light source to be saved in the current xml file 1460 XmlData.LIFCalib.Ray1Coord=ObjectData{1}.Coord; 1461 XmlData.LIFCalib.Ray2Coord=ObjectData{2}.Coord; 1462 XmlData.LIFCalib.RefLineCoord=ObjectData{3}.Coord; 1409 1463 1410 1464 %% display the current image in polar axes with origin at the illumination source … … 1414 1468 phys_polar=str2func('phys_polar'); 1415 1469 cd(currentdir) 1416 DataOut=phys_polar(UvData.Field,XmlData);1417 view_field(DataOut); 1418 1419 %% use the third line for reference luminosity 1470 XmlData.TransformInput.PolarCentre=[x0 y0]; 1471 DataPol=phys_polar(UvData.Field,XmlData);%transform the input image in polar coordinates with origin at the light source 1472 1473 %% use the third line for reference luminosity, renormalize the image intensity along each ray to get a uniform brightness along this line 1420 1474 if numel(val)==3 1421 1475 x_ref=linspace(ObjectData{3}.Coord(1,1),ObjectData{3}.Coord(2,1),10); … … 1424 1478 y_ref=y_ref-y0; 1425 1479 [theta_ref,r_ref] = cart2pol(x_ref,y_ref);%theta_ref and r_ref are the polar coordinates of the points on the line 1426 theta_ref=theta_ref*180/pi; 1427 figure 1480 theta_ref=theta_ref*180/pi;% theta_ref in radians 1481 figure(10) 1428 1482 plot(theta_ref,r_ref) 1429 azimuth_ima=linspace(DataOut.Coord_y(1),DataOut.Coord_y(2),size(DataOut.A,1));%profile of x index on the transformed image 1430 dist_source = interp1(theta_ref,r_ref,azimuth_ima); 1431 dist_source_pixel=round(size(DataOut.A,2)*(dist_source-DataOut.Coord_x(1))/(DataOut.Coord_x(2)-DataOut.Coord_x(1))); 1432 line_nan= isnan(dist_source_pixel); 1483 xlabel('azimuth(degrees)') 1484 ylabel('radius from light source') 1485 title('ref line in polar coordinates') 1486 azimuth_ima=linspace(DataPol.Coord_y(1),DataPol.Coord_y(2),size(DataPol.A,1));%array of angular index on the transformed image 1487 dist_source = interp1(theta_ref,r_ref,azimuth_ima);% get the polar position of the reference line 1488 dist_source_pixel=round(size(DataPol.A,2)*(dist_source-DataPol.Coord_x(1))/(DataPol.Coord_x(2)-DataPol.Coord_x(1))); 1489 line_nan= isnan(dist_source); 1433 1490 dist_source_pixel(line_nan)=1; 1434 1491 width=20; %number of pixels used for reference 1435 DataOut.A=double(DataOut.A); 1436 Anorm=zeros(size(DataOut.A)); 1437 Aval=mean(mean(DataOut.A)); 1438 for iline=1:size(DataOut.A,1) 1439 lum(iline)=mean(DataOut.A(iline,dist_source_pixel(iline):dist_source_pixel(iline)+width)); 1440 Anorm(iline,:)=uint16(Aval*DataOut.A(iline,:)/lum(iline)); 1492 DataPol.A=double(DataPol.A)-XmlData.LIFCalib.BlackOffset;% black background substracted 1493 Anorm=zeros(size(DataPol.A)); 1494 for iline=1:size(DataPol.A,1) 1495 lum(iline)=mean(DataPol.A(iline,dist_source_pixel(iline):dist_source_pixel(iline)+width));% average the luminosity on a band width lying on the reference line 1496 Anorm(iline,:)=DataPol.A(iline,:)/lum(iline);% for each ray (iline), renormalise the image by the brightness at the reference line 1441 1497 end 1442 1498 lum(line_nan)=NaN; 1443 figure 1444 plot(1:size(DataOut.A,1),lum) 1445 end 1446 ImaName=regexprep([get(handles.RootFile,'String') get(handles.FileIndex,'String')],'//',''); 1447 NewImageName=fullfile(get(handles.RootPath,'String'),'polar',[ImaName get(handles.FileExt,'String')]); 1448 imwrite(Anorm,NewImageName,'BitDepth',16) 1449 1450 %% record the origin in the xml file 1499 figure(11) 1500 plot(1:size(DataPol.A,1),lum) 1501 xlabel('azimuth (pixels)') 1502 ylabel('image brightness') 1503 %ImaName=regexprep([get(handles.RootFile,'String') get(handles.FileIndex,'String')],'//',''); 1504 NewImageName=fullfile(get(handles.RootPath,'String'),'LIF_ref.png'); 1505 imwrite(uint16(1000*Anorm),NewImageName,'BitDepth',16) 1506 DataPol.A=uint16(1000*Anorm); 1507 view_field(DataPol); 1508 end 1509 1510 % figure(12) 1511 % imagesc(DataPol.Coord_x,DataPol.Coord_y,AY,uint16(1000*Anorm)) 1512 %title('rescaled image') 1513 1514 %% record the calibration data in the xml file 1515 1451 1516 XmlFileName=find_imadoc(get(handles.RootPath,'String'),get(handles.SubDir,'String'),get(handles.RootFile,'String'),get(handles.FileExt,'String')); 1452 1517 answer=msgbox_uvmat('INPUT_Y-N','save the illumination origin in the current xml file?'); 1453 1518 if strcmp(answer,'Yes') 1454 1519 t=xmltree(XmlFileName); %read the file 1455 title =get(t,1,'name');1456 if ~strcmp(title ,'ImaDoc')1520 title_str=get(t,1,'name'); 1521 if ~strcmp(title_str,'ImaDoc') 1457 1522 msgbox_uvmat('ERROR','wrong xml file'); 1458 1523 return … … 1470 1535 return 1471 1536 end 1472 uid_illumination=find(t,'ImaDoc/ Illumination');1537 uid_illumination=find(t,'ImaDoc/LIFCalib'); 1473 1538 if isempty(uid_illumination) %if GeometryCalib does not already exists, create it 1474 [t,uid_illumination]=add(t,1,'element',' Illumination');1475 end 1476 uid_origin=find(t,'ImaDoc/ Illumination/Origin');1539 [t,uid_illumination]=add(t,1,'element','LIFCalib'); 1540 end 1541 uid_origin=find(t,'ImaDoc/LIFCalib/LightOrigin'); 1477 1542 if ~isempty(uid_origin) %if GeometryCalib does not already exists, create it 1478 1543 t=delete(t,uid_origin); 1479 1544 end 1480 % save the illumination origin 1481 t=struct2xml(XmlData.Illumination,t,uid_illumination); 1545 uid_line=find(t,'ImaDoc/LIFCalib/Ray1Coord'); 1546 if ~isempty(uid_line) %if GeometryCalib does not already exists, create it 1547 t=delete(t,uid_line); 1548 end 1549 uid_line=find(t,'ImaDoc/LIFCalib/Ray2Coord'); 1550 if ~isempty(uid_line) %if GeometryCalib does not already exists, create it 1551 t=delete(t,uid_line); 1552 end 1553 uid_line=find(t,'ImaDoc/LIFCalib/RefLineCoord'); 1554 if ~isempty(uid_line) %if GeometryCalib does not already exists, create it 1555 t=delete(t,uid_line); 1556 end 1557 uid_BlackOffset=find(t,'ImaDoc/LIFCalib/BlackOffset'); 1558 if ~isempty(uid_BlackOffset) %if GeometryCalib does not already exists, create it 1559 t=delete(t,uid_BlackOffset); 1560 end 1561 % save the LIF calibration data 1562 t=struct2xml(XmlData.LIFCalib,t,uid_illumination); 1482 1563 save(t,XmlFileName); 1483 1564 end 1565 1566 1567 %% get the image A*r 1568 [npy,npx]=size(DataPol.A); 1569 AX=DataPol.Coord_x; 1570 AY=DataPol.Coord_y; 1571 dX=(AX(2)-AX(1))/(npx-1); 1572 dY=(AY(1)-AY(2))/(npy-1);%mesh of new pixels 1573 r=AX(1)+[0:npx-1]*dX;%distance from laser 1574 rmat=ones(npy,1)*r; 1575 A=rmat.*Anorm; 1576 A=(A<=0).*ones(npy,npx)+(A>0).*A; %replaces zeros by ones 1577 A=log(A); % rA(r) should have a slope -gamma in x (radius from laser source) 1578 1579 %% find mask image: for use in uvmat 1580 % if ~isempty(DataIn.ZIndex) 1581 % num_level=DataIn.ZIndex;% needed to handle multilevel illumination 1582 % if ~exist('maskname','var') 1583 % huvmat=findobj(allchild(0),'tag','uvmat'); 1584 % hhuvmat=guidata(huvmat); 1585 % RootPath=get(hhuvmat.RootPath,'String'); 1586 % RootFile=get(hhuvmat.RootFile,'String'); 1587 % maskname=fullfile(RootPath,[RootFile '_23mask_' num2str(num_level) '.png']); 1588 % if test_1% PIV image used for correction 1589 % RootPath_1=get(hhuvmat.RootPath_1,'String'); 1590 % RootFile_1=get(hhuvmat.RootFile_1,'String'); 1591 % maskname_1=fullfile(RootPath_1,[RootFile_1 '_23mask_' num2str(num_level) '.png']); 1592 % end 1593 % end 1594 % Mask=DataIn; 1595 % if test_1 1596 % Mask_1=DataIn_1; 1597 % Mask.A=imread(maskname); 1598 % Mask_1.A=imread(maskname_1); 1599 % [Mask,Mask_1]=phys_polar(Mask,XmlData,Mask_1,XmlData_1); 1600 % Mask.A=Mask.A >= 200 & Mask_1.A >= 200 ; 1601 % % end 1602 % end 1603 % Mask.A=ones(size(DataPol.A)); 1604 [R,Y]=meshgrid(linspace(AX(1),AX(2),npx),linspace(AY(1),AY(2),npy)); 1605 [X,Y] = pol2cart(Y*pi/180,R); 1606 X=X+x0;% ny*nx matrix of phys cartesian coordiantes 1607 Y=Y+y0; 1608 gamma_coeff=zeros(1,npy); 1609 1610 NewImageName=fullfile(get(handles.RootPath,'String'),'LIF_ref_rescaled.png'); 1611 imwrite(uint16(1000*A),NewImageName,'BitDepth',16) 1612 return 1613 %loop on lines iY (angle in polar coordiantes) 1614 for iY=1:npy 1615 Theta_i=(AY(1)-(iY-1)*dY)*pi/180;%theta in radians 1616 %r_edge(iY)=(-100-XmlData.GeometryCalib.PolarCentre(1))/cos(Theta_i);%r_edge=radius of the edge (at X=-95) 1617 r_edge(iY) 1618 Aline=A(iY,:);% =r*image intensity 1619 Xline=X(iY,:); 1620 % adjust the fit region 1621 test_fit=X(iY,:)>10 & X(iY,:)<70 & A(iY,:)>0 & Mask.A(iY,:); 1622 ind_fit=find(test_fit);% indices of points used for the exponential fit 1623 if numel(ind_fit)<50 1624 Mask.A(iY,:)=0;% to exclude lines for which the line of rays is too short 1625 else 1626 p = polyfit(r(ind_fit),Aline(ind_fit),1); 1627 Aline_fit=polyval(p,r); 1628 if iY==round(npy/2) 1629 figure(12) 1630 plot(r(ind_fit),Aline(ind_fit),'b+',r,Aline_fit,'r') 1631 test_out=(Aline_fit-Aline)>0.5 1632 ind_out=find(test_out & test_fit) 1633 end 1634 test_out=(Aline_fit-Aline)>0.5;% eliminate points too dark compared to the fit, such that log(Afit/A)>0.5 1635 ind_out=find(test_out & test_fit); 1636 % if ~isempty(find(test_out & test_fit,1)) 1637 % Mask.A(iY,test_out)=0; 1638 if ~isempty(ind_out) 1639 Mask.A(iY,ind_out)=0; 1640 ind_fit=find(~test_out & test_fit); 1641 p = polyfit(r(ind_fit),Aline(ind_fit),1);%do new fit after elimination of the dark points 1642 end 1643 if -p(1) > 0.003 && -p(1) < 0.012 1644 gamma_coeff(iY)=-p(1); 1645 Aline_edge=polyval(p,r_edge(iY)); 1646 DataPol.A(iY,:)=exp(Aline_edge)/(r_edge(iY)); 1647 end 1648 1649 end 1650 end 1651 gamma_select=gamma_coeff(gamma_coeff~=0); 1652 gamma_select=filter(ones(1,20)/20,1,gamma_select);%smooth gamma_ceff over 20 pixels 1653 gamma_coeff(gamma_coeff~=0)=gamma_select; 1654 DataPol.A(Mask.A==0)=0; 1655 Ref.r_edge=r_edge; 1656 Ref.GammaCoeff=gamma_coeff 1657 1658 1659 1660 1661 1662 1663 1664 1665 1666 1484 1667 1485 1668 … … 2094 2277 end 2095 2278 end 2279 end 2280 if isfield(XmlDataRead, 'LIFCalib') 2281 XmlData.LIFCalib=XmlDataRead.LIFCalib; 2096 2282 end 2097 2283 end
Note: See TracChangeset
for help on using the changeset viewer.