%Read all notes throughout the code prior to running! %This code finds the yellow/green fluorescent boli in each of 8 cages after gavage %Go through JPG images in a folder and ID when cluster of pixels meets dyed poop criterion %Requires CAREFUL defining cages corner xy and exclusion of false positives %(incl cage corners) %Keep environment, light and camera settings constant ALWAYS! close all clear all %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Step one: Enter the pathname to your photos below %%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %pathname='C:\ENTER PATHNAME HERE' % Windows pathname='C:\JVVerhagen\PierceLab\Lab\Manuscripts\Gut motility odor\PICS\Week_2_Air\10_26_19_air' %pathname='/Users/kbaker/Documents/gut_transit/yellow_powder/Final_test/Week_1_Air/10_16_19_air' % MAC %threshold converting grey to black and white bw_th=0.45; %[0.45]=VALIDATED accurate and stable for fluoresc if based on green only; bw_th1=num2str(bw_th*100); %*100 avoids dot, avoiding file extension issues % min and max pixcount for boli; about 10-15 x 20-30 = ~200-450 pixels pixmin=100 %[100]=VALIDATED accurate and stable for fluoresc if based on green only; pixmax=600 %[600]=VALIDATED accurate and stable for fluoresc if based on green only; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Step two: Assign pixel values to the corners of the cages %%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Region Of Interest (ROI, the area defined by the 4 corners for each cage) %should NEVER overlap with top edge or corners of a non-black cage else FALSE POSITIVES %CAGES: %lower left cage = cage 1, upper left cage 2, lower right cage 7 %Thus, cage 1 is lower left corner; rest counted clockwise %CORNERS: %upper left corner first, upper right corner second, etc %Thus, first listed corner is upper leftcorner of the cage, rest listed clockwise %The pixel values can be determined by opening the image in ImageJ. %Outline the cage ROI, EDIT/SELECTION/PROPERTIES/LIST COORDINATES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % NOTE: SHIFT ALL ROIS with xshift and yshift %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% xshift= 0; yshift= 0; %cage1x=[x-value of upper left corner of cage 1; of upper right corner; of lower right corner; of lower left corner] %cage1y=[y-value of upper left corner of cage 1; of upper right corner; of lower right corner; of lower left corner] cage1x=[450+xshift; 1900+xshift; 1900+xshift; 450+xshift]; %start upper leftcorner of the cage, list clockwise cage1y=[2600+yshift; 2600+yshift; 3550+yshift; 3550+yshift]; cage2x=[200+xshift; 1000+xshift; 1000+xshift; 200+xshift]; %start upper leftcorner of the cage, list clockwise cage2y=[720+yshift; 720+yshift; 2340+yshift; 2340+yshift]; cage3x=[1200+xshift; 2300+xshift; 2300+xshift; 1200+xshift]; %start upper leftcorner of the cage, list clockwise cage3y=[720+yshift; 720+yshift; 2340+yshift; 2340+yshift]; cage4x=[2450+xshift; 3600+xshift; 3600+xshift; 2450+xshift]; %start upper leftcorner of the cage, list clockwise cage4y=[720+yshift; 720+yshift; 2340+yshift; 2340+yshift]; cage5x=[3700+xshift; 4900+xshift; 4900+xshift; 3700+xshift]; %start upper leftcorner of the cage, list clockwise cage5y=[720+yshift; 720+yshift; 2340+yshift; 2340+yshift]; cage6x=[5050+xshift; 5800+xshift; 5800+xshift; 5050+xshift]; %start upper leftcorner of the cage, list clockwise cage6y=[720+yshift; 720+yshift; 2340+yshift; 2340+yshift]; cage7x=[4150+xshift; 5700+xshift; 5700+xshift; 4150+xshift]; %start upper leftcorner of the cage, list clockwise cage7y=[2600+yshift; 2600+yshift; 3550+yshift; 3550+yshift]; cage8x=[2100+xshift; 3850+xshift; 3850+xshift; 2100+xshift]; %start upper leftcorner of the cage, list clockwise cage8y=[2600+yshift; 2600+yshift; 3550+yshift; 3550+yshift]; % from ColorThresholder app % Thresholds for fluorescent boli; % green & blue channel (emission) differentiates fluorescent boli vs rest Rrmin = 0; %[0] Grmin = 50; %[50] Brmin = 30; %[30] Rrmax = 50; %[50] Grmax = 255; %[255] Brmax = 255; %[255] %% creates a 2-D alpha shape of the points (x,y) cage1 = alphaShape(cage1x,cage1y); cage2 = alphaShape(cage2x,cage2y); cage3 = alphaShape(cage3x,cage3y); cage4 = alphaShape(cage4x,cage4y); cage5 = alphaShape(cage5x,cage5y); cage6 = alphaShape(cage6x,cage6y); cage7 = alphaShape(cage7x,cage7y); cage8 = alphaShape(cage8x,cage8y); %% FILES % rdpath = uigetdir(pathname, 'Choose the directory where Image files are located'); rdpath = pathname; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %NOTE: Ensure the correct rdpath1 is in use based on whether you are using a PC or a MAC. Either add or remove the % for the appropiate one. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% cd(rdpath) rdpath1=find(rdpath == '\');%FOR PC %rdpath1=find(rdpath == '/'); %FOR MAC foldername=[rdpath(rdpath1(end)+1:end)]; Names = dir('*.jpg'); fileSize = [Names(:).bytes]; selectedfiles = find(fileSize> 2000000); % & fileSize< 40000000); k=1; for k=1:length(selectedfiles) filename=Names(selectedfiles(k)).name; display(filename); I=imread(filename); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%NOTE:Figure(100) can be uncommented (remove "%") to see the the image that has been read into matlab %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %figure(100) %imagesc(I) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %NOTE: Uncomment figures 220,221,222 to visualize the red, green and %blue layers in the image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Ired=I(:,:,1); Igreen=I(:,:,2); Iblue=I(:,:,3); %figure(221);imshow(Igreen);title('GREEN LAYER ONLY') %figure(220);imshow(Ired);title('RED LAYER ONLY') %figure(222);imshow(Iblue); title('BLUE LAYER ONLY') %% remove mice & noise by identifying continious regions of certain size clear CC numPixels idx idx2 idx3 BWI %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%NOTE:Below can be uncommented to see the black and white image % BWI1 = rgb2gray(I); % figure(1001);imshow(BWI1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %BWI = im2bw(rgb2gray(I), bw_th); %BWI = im2bw(Iblue, bw_th); %size based on blue only BWI = im2bw(Igreen, bw_th); %size based on green only %figure(152);imshow(BWI); %BW only above high greenness BWIc = imcomplement(BWI); BWIc = BWI; CC = bwconncomp(BWIc); %list of fully connected pixels, % k-th element in the cell array is a vector containing the linear indices of the pixels in the k-th object. numPixels = cellfun(@numel,CC.PixelIdxList); %size of each connected element %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%NOTE:Below can be uncommented to see the black and white image %figure(1012);hist(numPixels); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% idx = find(numPixels>pixmax); %mice etc; index in CC of connected elements w that size idx2 = find(numPixelspixmin & numPixels0); %boli indices if InCage_ind >0 clear id4 inds4 boli_r boli_g boli_b boli_rgb_ranged for id4=1:length(InCage_ind) inds4=CC.PixelIdxList{idx3(InCage_ind(id4))}; %pixels per boli boli_r(id4)=mean(Ired(inds4)); %red value of all pixels in a boli boli_g(id4)=mean(Igreen(inds4)); boli_b(id4)=mean(Iblue(inds4)); %Boli within RGB range; separate mins and max boli_rgb_ranged(id4) = (boli_r(id4) >= Rrmin) & (boli_r(id4) <= Rrmax) & ... (boli_g(id4) >= Grmin) & (boli_g(id4) <= Grmax) & ... (boli_b(id4) >= Brmin) & (boli_b(id4) <= Brmax) ; end Boli_rgb_total(k,numcage)=sum(boli_rgb_ranged); else %Boli2_total(k, numcage)=0; Boli_rgb_total(k,numcage)=0; end end %%%%%%%%%%%%% show per jpg for troubleshooting %figure(302); plot(boli_sums(k,:));hold on; plot(Boli_rgb_total(k,:)); hold off %xlabel('Cages'); ylabel('Num Poops'); %legend('by size and green', 'boli also by color') %Summary(k+1,[1:5])={[Names(selectedfiles(k)).name boli_sums Boli_rgb_total]}; end %% PLOTTING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %NOTE: %Figure 300 with give a image of the pixel outline of the cages and is saved %Figure 311 is a graph to show the number of boli change over time and it is saved %Figure 320 is the origional cage image overlayed with cage regions %Figure 321 is the origional cage image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure(300); h=figure(300); axis on; %needed for final subplots plot(cage1);axis off; hold on; plot(cage2);axis off; hold on; plot(cage3);axis off; hold on; plot(cage4);axis off; hold on; plot(cage5);axis off; hold on; plot(cage6);axis off; hold on; plot(cage7);axis off; hold on; plot(cage8);axis on; xlim([0 size(I,2)]); ylim([0 size(I,1)]); hold off saveas(h,'cages','jpg') for l=1:8 figure(311); subplot(3,3,l) plot(boli_sums(:,l),'.');hold on; plot(Boli_rgb_total(:,l),'+'); hold off xlabel('Time (img)'); ylabel('Num Boli'); h11=legend('by size', 'by size & color'); set(h11, 'Position', [0.05, 0.9, 0.25, 0.1]) cagenum=num2str(l); title(['Cage ' cagenum]) end %%%%%%%%%%%%% for troubleshooting % subplot(3,3,8); imagesc(I);imagenum=num2str(k); title(['Time ' imagenum 'min']) % show overlayed images as a subplot subplot(3,3,9); imagesc(I), hold on cageim=imread('cages.jpg'); cageim1=flipdim(cageim,1); cageim2= imresize(cageim1,[size(I,1) size(I,2)]); %hImg = imshow(cageim2); set(hImg, 'AlphaData', 0.6) %get corners for cutout cageim3=cageim2(690:3450, 780:5440); cageim4= imresize(cageim3,[size(I,1) size(I,2)]); hImg = imshow(cageim4); set(hImg, 'AlphaData', 0.6) %Saves figure 311 to folder h2=figure(311); set(gcf,'units','normalized','outerposition',[0 0 1 1]); pixmin1=num2str(pixmin);pixmax1=num2str(pixmax); saveas(h2,['gavage ' bw_th1 ' ' pixmin1 '-' pixmax1 ' ' foldername],'jpg'); %overlay of cages with cage region figure(320) imagesc(I), hold on hImg = imshow(cageim4); set(hImg, 'AlphaData', 0.6) %origional image figure(321) imagesc(I) clear I cageim2 CC cageim3 cageim4 Iblue Igreen Ired %calculate the time of the first boli clear firstboli firstbolus o for o=1:size(Boli_rgb_total,2) firstboli=min(find(Boli_rgb_total(:,o)>0)); if isempty(firstboli) firstboli=100; %dummy, no boli in this cage end firstbolus(o)=firstboli; end open Boli_rgb_total open firstbolus %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%For PC - comment out accordingly %%%%%%%%%%%%%%%%%%%%%%%%%%%%% filenamexls = ['allboli_rgb_total ' bw_th1 ' ' pixmin1 '-' pixmax1 ' ' foldername '.xlsx']; writematrix(Boli_rgb_total,filenamexls) filenamexls = ['firstbolus ' bw_th1 ' ' pixmin1 '-' pixmax1 ' ' foldername '.xlsx']; writematrix(firstbolus,filenamexls) %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%For Mac - comment out accordingly %%%%%%%%%%%%%%%%%%%%%%%%%%%%% dlmwrite('allboli.xls', Boli_rgb_total,'\t') dlmwrite('firstbolus.xls', firstbolus,'\t') %% CURVE FITS clear session_dynamics fit for cagei= 1:size(Poops_rgb_total,2) [xData, yData1] = prepareCurveData( 1:size(Poops_rgb_total,1), Poops_rgb_total(:,cagei) ); yData=smooth(yData1,5); ft = fittype('rise_amplitude ./ (exp(-4.4*(x - onset_time) ./ rise_time) + 1)', 'independent', 'x', 'dependent', 'y' ); opts = fitoptions( 'Method', 'NonlinearLeastSquares','Lower', [0 0 0],'Upper', [80 20 80] ); [fitresult, gof] = fit(xData, yData, ft, opts); figure(391) plot(fitresult,xData, yData) coeffnames(fitresult) [coeffsfit]= coeffvalues(fitresult); onset_tfit(cagei)=coeffsfit(1); rise_ampfit(cagei)=coeffsfit(2); rise_tfit(cagei)=coeffsfit(3); rsqsfit(cagei)= gof.rsquare; end session_dynamics=[onset_tfit; rise_ampfit; rise_tfit; rsqsfit] open session_dynamics % Excluding empty cage 6 all 3 groups, empty cage 2 group 1: 20 mice of Fig2C % Day 1 Day 2 % avg sd avg sd % first bolus time 263.5 140.9 222.5 81.4 % onset time 254.6 109.3 227.3 88.3 % rise time 116.8 122.0 69.4 82.8 % number boli 2.4 3.2 2.2 1.7 % fit rsq 0.87 0.20 0.84 0.25 % The coefficient of variation (sd/avg) was reduced by 19 on day 1 (0.53 vs 0.43), % but increased by 8% on day 2 (0.36 vs 0.39). % This approach also yields new insights in dynamics via the rise times, being 117 and 69 min, respectively. % The mean number of boli were 2.4 and 2.2. % Quality of the fits was high fitting 84-85 of the variance of the data. %% VALIDATION % % validation=csvread(validationfile); %bolus counts in 8 cages (8 collumns) x images (rows) % validation_vect=reshape(validation,[],[1]); % % Poops_rgb_total_vect=reshape(Poops_rgb_total,[],[1]); % Numcagepics=length(Poops_rgb_total_vect); % valid_vect_r=corr(validation_vect, Poops_rgb_total_vect); % valid_vect_diff_m=mean(validation_vect-Poops_rgb_total_vect); % valid_vect_diff_m_abs=mean(abs(validation_vect-Poops_rgb_total_vect)); % valid_vect_diff_s=std(validation_vect-Poops_rgb_total_vect); % % display('Validation stats: n; correl; mean diff; mean abs diff; sd diff') % valid_stats=[Numcagepics valid_vect_r valid_vect_diff_m valid_vect_diff_m_abs valid_vect_diff_s] % % figure(342) % %scatter(validation_vect, Poops_rgb_total_vect) % plot(Poops_rgb_total(:, 1),'b'); hold on; % plot(validation(:,1),'b+');hold on; % plot(Poops_rgb_total(:, 2),'k'); hold on; % plot(validation(:,2),'k+');hold on; % plot(Poops_rgb_total(:, 3),'r'); hold on; % plot(validation(:,3),'r+');hold on; % plot(Poops_rgb_total(:, 4),'g'); hold on; % plot(validation(:,4),'g+');hold on; % plot(Poops_rgb_total(:, 5),'y'); hold on; % plot(validation(:,5),'y+');hold on; % plot(Poops_rgb_total(:, 6),'w'); hold on; % plot(validation(:,6),'w+');hold on; % plot(Poops_rgb_total(:, 7),'m'); hold on; % plot(validation(:,7),'m+');hold on; % plot(Poops_rgb_total(:, 8),'c'); hold on; % plot(validation(:,8),'c+');hold on; % hold off % % h3=figure(342); % xlabel('Time (img)'); ylabel('Num Poops'); % set(gcf,'units','normalized','outerposition',[0 0 1 1]); % % valid_numcagepics1=num2str(Numcagepics); % valid_vect_r1=num2str(valid_vect_r); % valid_vect_diff_m1=num2str(valid_vect_diff_m); % valid_vect_diff_m_abs1=num2str(valid_vect_diff_m_abs); % valid_vect_diff_s1=num2str(valid_vect_diff_s); % title(['n=' valid_numcagepics1 '; Validation correl: ' valid_vect_r1 '; mean diff: ' valid_vect_diff_m1 '; mean abs diff: ' valid_vect_diff_m_abs1 '; SD: ' valid_vect_diff_s1 ]) % legend('mouse 1 auto','valid','mouse 2 auto', 'valid','mouse 3 auto','valid', 'mouse 4 auto','valid','mouse 5 auto','valid','mouse 6 auto','valid','mouse 7 auto','valid','mouse 8 auto','valid') % saveas(h3,['Validation ' bw_th1 ' ' pixmin1 '-' pixmax1 ' ' foldername],'jpg');