Changeset 1049
 Timestamp:
 Jun 5, 2018, 10:13:08 PM (5 years ago)
 Location:
 trunk/src
 Files:

 3 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/find_field_cells.m
r1048 r1049 475 475 for ivar=ind_coord_y 476 476 DimCell=Data.VarDimName{ivar}; 477 if numel(DimCell )==1 && strcmp(DimCell_x{1},DimCell{1})477 if numel(DimCell_x)==1 && strcmp(DimCell_x{1},DimCell{1}) 478 478 y_nbre(icell)=y_nbre(icell)+1; 479 479 Cell1DPlot{icell}.YIndex(y_nbre(icell))=ivar; 
trunk/src/series/extract_multitif_special.m
r1048 r1049 144 144 % end 145 145 %ImagesPerLevel=size(XmlInput.Time,2)1;%100;%use the xmlinformation to get the nbre of j indices 146 Dtj=0.05;% time interval between frames 146 147 ImagesPerLevel=450;% total number of images per position, ImagesPerLevelNbj images skiiped during motion between two positions 147 148 Nbj=400; %Nbre of images kept at a given position 149 Dti=Dtj*ImagesPerLevel; 150 NbLevel=11; 151 NbScan=3; 152 TimeReturn=20; %time needed to return back to the first position (in sec) 153 NbSkippedReturn=round(TimeReturn/Dtj); 148 154 %% create the xml file of PCO camera if it does not exist 149 155 Newxml=fullfile(Param.InputTable{1,1},[Param.InputTable{1,2} '.xml']); … … 154 160 % XmlInput.Camera.BurstTiming.Time=1;% for 180 155 161 XmlInput.Camera.BurstTiming.Time=0;% for 200 156 XmlInput.Camera.BurstTiming.Dtj=0.05; 157 XmlInput.Camera.BurstTiming.NbDtj=199; 158 XmlInput.Camera.BurstTiming.Dti=12.5; 159 XmlInput.Camera.BurstTiming.NbDti=10; 160 XmlInput.Camera.BurstTiming.Dtk=157; 161 XmlInput.Camera.BurstTiming.NbDtk=2; 162 XmlInput.Camera.BurstTiming.Dtj=Dtj; 163 XmlInput.Camera.BurstTiming.NbDtj=Nbj1; 164 XmlInput.Camera.BurstTiming.Dti=Dti; 165 XmlInput.Camera.BurstTiming.NbDti=NbLevel1; 166 XmlInput.Camera.BurstTiming.Dtk=Dti*Nb 167 Level+TimeReturn; 168 XmlInput.Camera.BurstTiming.NbDtk=NbScan1; 162 169 %XmlInput=rmfield(XmlInput,'Time'); 163 170 %XmlInput=rmfield(XmlInput,'TimeUnit'); … … 198 205 i_index=fix((count1)/ImagesPerLevel)+1; 199 206 j_index=mod(count1,ImagesPerLevel)+1; 200 elseif count<=ImagesPerLevel* 11+400 %skip 400 images during return207 elseif count<=ImagesPerLevel*NbLevel+400 %skip 400 images during return 201 208 checkkeep=0; 202 elseif count<=ImagesPerLevel*2 2+400 % =5930 second scan209 elseif count<=ImagesPerLevel*2*NbLevel+400 % =5930 second scan 203 210 i_index=fix((count401)/ImagesPerLevel)+1; 204 211 j_index=mod(count401,ImagesPerLevel)+1; 205 elseif count<=ImagesPerLevel*2 2+800 %skip images during second return, from 2763 to 3167212 elseif count<=ImagesPerLevel*2*NbLevel+800 %skip images during second return, from 2763 to 3167 206 213 checkkeep=0; 207 elseif count<=ImagesPerLevel*3 3+800 % =5930 third scan214 elseif count<=ImagesPerLevel*3*NbLevel+800 % =5930 third scan 208 215 i_index=fix((count801)/ImagesPerLevel)+1; 209 216 j_index=mod(count801,ImagesPerLevel)+1; 
trunk/src/uvmat.m
r1048 r1049 1371 1371 1372 1372 %% read the lines previously used for LIF calibration if available, and display them 1373 ListObjectName=get(handles.ListObject,'String'); 1373 1374 if isfield(XmlData,'LIFCalib') 1374 1375 if isfield(XmlData.LIFCalib,'Ray1Coord') && isfield(XmlData.LIFCalib,'Ray2Coord') && isfield(XmlData.LIFCalib,'RefLineCoord') 1376 ListObjectName(1:4)={'plane';'Ray1';'Ray2';'RefLine'}; 1375 1377 data.Name='Ray1'; 1376 1378 data.Type='line'; … … 1388 1390 UvData.ProjObject{4}.DisplayHandle.uvmat=plot_object(UvData.ProjObject{4},UvData.ProjObject{1},handles.PlotAxes,'b'); 1389 1391 end 1390 if isfield(XmlData.LIFCalib,'MaskPolygon') 1392 if isfield(XmlData.LIFCalib,'MaskPolygonCoord') 1393 ListObjectName{5}='MaskPolygon'; 1391 1394 data.Name='MaskPolygon'; 1392 1395 data.Coord=XmlData.LIFCalib.MaskPolygonCoord; 1393 1396 UvData.ProjObject{5}=data; 1394 1397 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 end 1399 end 1400 set(handles.ListObject,'String',ListObjectName) 1398 1401 set(handles.ListObject,'Value',1) 1399 1402 set(handles.uvmat,'UserData',UvData); … … 1401 1404 %% read lines currently drawn 1402 1405 ListObj=UvData.ProjObject;% get the current list of projection objects 1403 select =zeros(1,numel(ListObj));1406 select_line=zeros(1,numel(ListObj)); 1404 1407 select_mask=zeros(1,numel(ListObj)); 1405 1408 for iobj=1:numel(ListObj); 1406 1409 if isfield(ListObj{iobj},'Type') && strcmp(ListObj{iobj}.Type,'line') 1407 select (iobj)=1;% select the lines among the projection objects1410 select_line(iobj)=1;% select the lines among the projection objects 1408 1411 end 1409 1412 if isfield(ListObj{iobj},'Type') && strcmp(ListObj{iobj}.Type,'polygon') … … 1411 1414 end 1412 1415 end 1413 val=find(select); 1414 if numel(val)<2 1416 if numel(find(select_line))<2 1415 1417 msgbox_uvmat('ERROR',{'light rays must be defined by at least two lines created by Projection object/line in the menu bar'; ... 1416 1418 'use a third line to get a reference luminosity profile accross the illumination beam'}); 1417 1419 return 1418 1420 else 1419 ObjectData=UvData.ProjObject(val); 1420 for iobj=1:length(ObjectData) 1421 xA(iobj)=ObjectData{iobj}.Coord(1,1); 1422 yA(iobj)=ObjectData{iobj}.Coord(1,2); 1423 xB(iobj)=ObjectData{iobj}.Coord(2,1); 1424 yB(iobj)=ObjectData{iobj}.Coord(2,2); 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 1421 LineData=UvData.ProjObject(find(select_line)); 1422 for iobj=1:length(LineData) 1423 xA(iobj)=LineData{iobj}.Coord(1,1); 1424 yA(iobj)=LineData{iobj}.Coord(1,2); 1425 xB(iobj)=LineData{iobj}.Coord(2,1); 1426 yB(iobj)=LineData{iobj}.Coord(2,2); 1427 end 1428 end 1429 1434 1430 1435 1431 %% set the image offset … … 1442 1438 return 1443 1439 else 1444 XmlData.LIFCalib.BlackOffset= answer;% image value for black background, to be determined by taking images with a cover on the objective lens1440 XmlData.LIFCalib.BlackOffset=str2num(answer) ;% image value for black background, to be determined by taking images with a cover on the objective lens 1445 1441 end 1446 1442 … … 1458 1454 y0=((y3y4)*(x1*y2y1*x2)(y1y2)*(x3*y4y3*x4))/D; 1459 1455 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;1456 XmlData.LIFCalib.Ray1Coord=LineData{1}.Coord; 1457 XmlData.LIFCalib.Ray2Coord=LineData{2}.Coord; 1458 XmlData.LIFCalib.RefLineCoord=LineData{3}.Coord; 1463 1459 1464 1460 %% display the current image in polar axes with origin at the illumination source … … 1472 1468 1473 1469 %% use the third line for reference luminosity, renormalize the image intensity along each ray to get a uniform brightness along this line 1474 if numel( val)==31475 x_ref=linspace( ObjectData{3}.Coord(1,1),ObjectData{3}.Coord(2,1),10);1476 y_ref=linspace( ObjectData{3}.Coord(1,2),ObjectData{3}.Coord(2,2),10);1470 if numel(find(select_line))==3 1471 x_ref=linspace(LineData{3}.Coord(1,1),LineData{3}.Coord(2,1),10); 1472 y_ref=linspace(LineData{3}.Coord(1,2),LineData{3}.Coord(2,2),10); 1477 1473 x_ref=x_refx0; 1478 1474 y_ref=y_refy0; … … 1501 1497 xlabel('azimuth (pixels)') 1502 1498 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); 1499 title('ref brightness profile') 1500 end 1501 1502 %% get the image A*r 1503 [npy,npx]=size(Anorm); 1504 AX=DataPol.Coord_x; 1505 AY=DataPol.Coord_y; 1506 dX=(AX(2)AX(1))/(npx1); 1507 dY=(AY(1)AY(2))/(npy1);%mesh of new pixels 1508 [R,Theta]=meshgrid(linspace(AX(1),AX(2),npx),linspace(AY(1),AY(2),npy)); 1509 A=R.*Anorm; 1510 A=(A<=0).*ones(npy,npx)+(A>0).*A; %replaces zeros by ones 1511 A=log(A); % rA(r) should have a slope gamma in x (radius from laser source) 1512 1513 %% get the exponential decay law 1514 index_mask=find(select_mask); 1515 if ~isempty(index_mask) 1516 if numel(index_mask)>1 1517 msgbox_uvmat('ERROR',{'only one mask polygon accepted'}); 1518 return 1519 else 1520 MaskData=UvData.ProjObject{index_mask}; 1521 XmlData.LIFCalib.MaskPolygonCoord=MaskData.Coord; 1522 end 1523 end 1524 1525 %loop on lines iY (angle in polar coordiantes) 1526 gamma_coeff=zeros(1,npy); 1527 for iY=1:npy 1528 ALine=A(iY,:); 1529 RLine=R(iY,:); 1530 ThetaLine=Theta(iY,:)*pi/180; 1531 [XLine,YLine] = pol2cart(ThetaLine,RLine); 1532 XLine=XLine+x0; 1533 YLine=YLine+y0; 1534 if ~isempty(index_mask) 1535 ind_good=inpolygon(XLine,YLine,MaskData.Coord(:,1),MaskData.Coord(:,2)); 1536 if numel(find(ind_good))>100 1537 ALine=ALine(ind_good); 1538 RLine=RLine(ind_good); 1539 else 1540 continue 1541 end 1542 end 1543 p = polyfit(RLine,ALine,1); 1544 gamma_coeff(iY)=p(1); 1545 end 1546 1547 %% show and store the rescaled image 1548 NewImageName=fullfile(get(handles.RootPath,'String'),'LIF_ref_rescaled.png'); 1549 DataPol.A=uint16(1000*A); 1507 1550 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')1551 imwrite(DataPol.A,NewImageName,'BitDepth',16) 1552 figure(13) 1553 plot(Theta(:,1),100*gamma_coeff) 1554 xlabel('angle(degree)') 1555 ylabel('decay coeff(m1)') 1513 1556 1514 1557 %% record the calibration data in the xml file 1515 1516 1558 XmlFileName=find_imadoc(get(handles.RootPath,'String'),get(handles.SubDir,'String'),get(handles.RootFile,'String'),get(handles.FileExt,'String')); 1517 1559 answer=msgbox_uvmat('INPUT_YN','save the illumination origin in the current xml file?'); … … 1551 1593 t=delete(t,uid_line); 1552 1594 end 1553 uid_ line=find(t,'ImaDoc/LIFCalib/RefLineCoord');1554 if ~isempty(uid_ line) %if GeometryCalib does not already exists, create it1555 t=delete(t,uid_ line);1595 uid_mask=find(t,'ImaDoc/LIFCalib/MaskPolygonCoord'); 1596 if ~isempty(uid_mask) %if GeometryCalib does not already exists, create it 1597 t=delete(t,uid_mask); 1556 1598 end 1557 1599 uid_BlackOffset=find(t,'ImaDoc/LIFCalib/BlackOffset'); … … 1563 1605 save(t,XmlFileName); 1564 1606 end 1565 1566 1567 %% get the image A*r1568 [npy,npx]=size(DataPol.A);1569 AX=DataPol.Coord_x;1570 AY=DataPol.Coord_y;1571 dX=(AX(2)AX(1))/(npx1);1572 dY=(AY(1)AY(2))/(npy1);%mesh of new pixels1573 r=AX(1)+[0:npx1]*dX;%distance from laser1574 rmat=ones(npy,1)*r;1575 A=rmat.*Anorm;1576 A=(A<=0).*ones(npy,npx)+(A>0).*A; %replaces zeros by ones1577 A=log(A); % rA(r) should have a slope gamma in x (radius from laser source)1578 1579 %% find mask image: for use in uvmat1580 % if ~isempty(DataIn.ZIndex)1581 % num_level=DataIn.ZIndex;% needed to handle multilevel illumination1582 % 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 correction1589 % 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 % end1593 % end1594 % Mask=DataIn;1595 % if test_11596 % 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 % % end1602 % end1603 % 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 coordiantes1607 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 return1613 %loop on lines iY (angle in polar coordiantes)1614 for iY=1:npy1615 Theta_i=(AY(1)(iY1)*dY)*pi/180;%theta in radians1616 %r_edge(iY)=(100XmlData.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 intensity1619 Xline=X(iY,:);1620 % adjust the fit region1621 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 fit1623 if numel(ind_fit)<501624 Mask.A(iY,:)=0;% to exclude lines for which the line of rays is too short1625 else1626 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_fitAline)>0.51632 ind_out=find(test_out & test_fit)1633 end1634 test_out=(Aline_fitAline)>0.5;% eliminate points too dark compared to the fit, such that log(Afit/A)>0.51635 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 points1642 end1643 if p(1) > 0.003 && p(1) < 0.0121644 gamma_coeff(iY)=p(1);1645 Aline_edge=polyval(p,r_edge(iY));1646 DataPol.A(iY,:)=exp(Aline_edge)/(r_edge(iY));1647 end1648 1649 end1650 end1651 gamma_select=gamma_coeff(gamma_coeff~=0);1652 gamma_select=filter(ones(1,20)/20,1,gamma_select);%smooth gamma_ceff over 20 pixels1653 gamma_coeff(gamma_coeff~=0)=gamma_select;1654 DataPol.A(Mask.A==0)=0;1655 Ref.r_edge=r_edge;1656 Ref.GammaCoeff=gamma_coeff1657 1658 1659 1660 1661 1662 1663 1664 1607 1665 1608
Note: See TracChangeset
for help on using the changeset viewer.