diff --git a/CavBubbleDiameter.m b/CavBubbleDiameter.m
new file mode 100644
index 0000000..edcaa9f
--- /dev/null
+++ b/CavBubbleDiameter.m
@@ -0,0 +1,85 @@
+clc;
+clear;
+
+basedir='...\Visualization\';
+
+LP=[300 400 500 600 700 800 900];
+magf=35;% (pixels/mm)
+
+LPmJ=0.1122*LP-20.61;
+
+Nl=7;
+Nr=5;
+
+x0=zeros(Nl,Nr);
+y0=zeros(Nl,Nr);
+bubble_diameter=zeros(Nl,Nr);
+bubble_area=zeros(Nl,Nr);
+cropy=200:640;
+
+
+for i=1:Nl
+    
+    for j=1:Nr
+        
+        
+        % Read frame
+        A=imread([basedir,'LP',num2str(LP(i),'%03.0f'),filesep,'run',int2str(j),'\cam1__C001H001S0001000010.tif']);
+            
+        img=A(end:-1:1,:);
+        
+        
+        % crop image
+        img = img(cropy,:);
+        % flip the image upside down
+        img = flipud(img);
+        
+        img(img>100)=255;
+%           % converts image to black and white
+          BW = imbinarize(img,0.99);
+        
+        
+        % only pixels that are connected to another # pixels are kept
+        BWred = bwareaopen(~BW,100);
+        % fills in holes of the image
+        BWfill=imfill(BWred,'holes');
+        
+        
+        
+        % finds the boundaries of the image
+        [B,labelmat,N,~]= bwboundaries(BWfill);
+        [~,ind]=max(cellfun('length',B));
+        bound = B{ind};
+        [N ind]
+        % calculate geometrical properties of the bubble
+        indx=find(labelmat==ind);
+        bubble_area(i,j)=length(indx);      
+        bubble_diameter(i,j) = sqrt(bubble_area(i,j) * 4/pi); % * pixel_pitch/magnification;
+        
+        if bubble_diameter(i,j)<=30
+            ind=ind+1;
+            bound = B{ind};
+            indx=find(labelmat==ind);
+            bubble_area(i,j)=length(indx);      
+            bubble_diameter(i,j) = sqrt(bubble_area(i,j) * 4/pi); % * pixel_pitch/magnification;
+        end
+        
+        btemp=bound;
+        % calculate centroid
+        x0(i,j) = mean(btemp(:,2));
+        y0(i,j) = mean(btemp(:,1));
+        
+        
+        
+%         [x0 y0 bubble_diameter]
+%%
+%         figure(1);imshow(BWfill(:,:));hold on;
+% %         btemp=bound{1};
+%         plot(btemp(:,2),btemp(:,1),'r-','LineWidth',2);
+%         hold off;
+%         pause(0.3);
+
+    end
+end
+
+