# Find the Length of a Pendulum in Motion

This example shows you how to calculate the length of a pendulum in motion by segmenting the pendulum and calculating the center of the pendulum across many image frames, then fitting the center coordinates to the equation of a circle.

Load the pendulum data from a MAT file. The file contains two variables. The `frames` variable contains the video as a 4-D numeric array. The `rect` variable contains a region of interest that bounds the extent of the pendulum.

```load pendulum; ```

Preview the video using `implay`. You can see that the pendulum is swinging in the upper half of each frame.

```implay(frames); ``` ### Step 2: Select Region where Pendulum is Swinging

To improve the efficiency of the segmentation, create a new series of frames that contains only the region where the pendulum is swinging.

First perform `imcrop` on one frame, specifyiing the cropping region as the preloaded `rect` variable. Next, initialize an array to store all of the cropped frames, based on the size of the first cropped frame. Finally, loop through all of the frames, crop each one using the same cropping region, and save the result to the initialized array.

```first_frame = frames(:,:,:,1); first_region = imcrop(first_frame,rect); nFrames = size(frames,4); frame_regions = repmat(uint8(0), [size(first_region) nFrames]); for count = 1:nFrames frame_regions(:,:,:,count) = imcrop(frames(:,:,:,count),rect); end imshow(frame_regions(:,:,:,1)) ``` ### Step 3: Segment the Pendulum in Each Frame

Notice that the pendulum is much darker than the background. You can segment the pendulum in each frame by converting the frame to grayscale, thresholding it using `imbinarize`, and removing background structures using `imopen` and `imclearborder`.

Initialize a binary array to contain the segmented pendulum frames.

```seg_pend = false([size(first_region,1) size(first_region,2) nFrames]); centroids = zeros(nFrames,2); se_disk = strel("disk",3); ```

Loop through all of the frames and perform the segmentation and postprocessing. In a montage, display the original frame and the segmentation result.

```for count = 1:nFrames fr = frame_regions(:,:,:,count); gfr = im2gray(fr); gfr = imcomplement(gfr); bw = imbinarize(gfr,0.7); % Threshold is determined experimentally bw = imopen(bw,se_disk); bw = imclearborder(bw); seg_pend(:,:,count) = bw; montage({fr,labeloverlay(gfr,bw)}); pause(0.2) end ``` ### Step 4: Find the Center of the Segmented Pendulum in Each Frame

You can see that the shape of the pendulum varied in different frames. This is not a serious issue because you just need its center. You will use the pendulum centers to find the length of the pendulum.

Use `regionprops` to calculate the center of the pendulum.

```pend_centers = zeros(nFrames,2); for count = 1:nFrames property = regionprops(seg_pend(:,:,count),"Centroid"); pend_centers(count,:) = property.Centroid; end ```

Display the pendulum centers using `plot`.

```x = pend_centers(:,1); y = pend_centers(:,2); figure plot(x,y,"m.") axis ij axis equal hold on; xlabel("x"); ylabel("y"); title("Pendulum Centers"); ``` ### Step 5: Calculate Radius by Fitting a Circle Through Pendulum Centers

Rewrite the basic equation of a circle:

`(x-xc)^2 + (y-yc)^2 = radius^2`

where `(xc,yc)` is the center, in terms of parameters `a`, `b`, `c` as

`x^2 + y^2 + a*x + b*y + c = 0`

where `a = -2*xc`, `b = -2*yc`, and `c = xc^2 + yc^2 - radius^2`.

You can solve for parameters `a`, `b`, and `c` using the least squares method. Rewrite the above equation as

`a*x + b*y + c = -(x^2 + y^2)`

which can also be rewritten as

`[x y 1] * [a;b;c] = -x^2 - y^2`.

Solve this equation using the backslash(`\`) operator.

The circle radius is the length of the pendulum in pixels.

```abc = [x y ones(length(x),1)] \ -(x.^2 + y.^2); a = abc(1); b = abc(2); c = abc(3); xc = -a/2; yc = -b/2; circle_radius = sqrt((xc^2 + yc^2) - c); pendulum_length = round(circle_radius) ```
```pendulum_length = 253 ```

Superimpose the circle and circle center on the plot of pendulum centers.

```circle_theta = pi/3:0.01:pi*2/3; x_fit = circle_radius*cos(circle_theta)+xc; y_fit = circle_radius*sin(circle_theta)+yc; plot(x_fit,y_fit,"b-"); plot(xc,yc,"bx","LineWidth",2); plot([xc x(1)],[yc y(1)],"b-"); text(xc-110,yc+100,sprintf("Pendulum length = %d pixels",pendulum_length)); ``` 