PROGRAM DETAILS

The entire software is written in python. The source code and the documentation (html format using sphinx) can be found here (https://github.com/kanishk-ux/Tethered-cell-analysis.git)
It uses various functions defined in an open source computer vision and machine learning software library called OpenCV(about OpenCV https://opencv.org/about/). Other python libraries used are NumPy, sklearn (linear_model), matplotlib, scipy to make various calculations as per the requirements. (Details about these can be found here https://pythontips.com/2018/06/03/top-14-most-famous-python-libraries-frameworks/).
 
Some features of the program : 
Has a user friendly GUI  to choose the folder to analyse and give initial data for the program to run.
The program is very fast (Can analyse 500 frames in less than 20 seconds after the cell to be analysed is selected.)
The program shows a preview for the user so that it is easy for him/her to select the desired cell for analysis since a single video can have multiple rotating tethered cells.
The program can be slowed down by a user defined delay value so that the user can see if the displayed regression fit for each frame of the cell is correct or not. 
Calculates the following values using the raw image stacks generated by imagej after microscopy:
1. Angular velocity (and the mean and standard deviation)
2. Frequency, Average CCW. CW frequency
3. Number of CCW, CW rotations and total time spent in each direction. ( also the number of    CCW to CW switches and CW to CCW switches )
4. Length of the cell, Trace diameter.
5. Center of rotation of the cell and distance between the COM and center of rotation.
6. Cumulative angle.
7. Plots the COM data for each frame, Cumulative angle vs Frame count, Frequency vs time, fitting circles to the data using three methods and calculating residual values and choosing the best method. Dynamically adjusting for the center of rotation.
8. Saving all plots and excel files generated from the images in a single folder at a desired path.

What is the thresholding algorithm used in the program? (different filters and how they are different)
 
There are many forms of image segmentation like Clustering, Compression, Edge detection, Region-growing, Graph partitioning, watershed and the most basic type of image segmentation: thresholding.
The function used is cv2.threshold. First argument is the source image, which should be a grayscale image. Second argument is the threshold value which is used to classify the pixel values. Third argument is the maxVal which represents the value to be given if pixel value is more than (sometimes less than) the threshold value. OpenCV provides different styles of thresholding and it is decided by the fourth parameter of the function. Here we use cv2.THRESH_BINARY.
All thresholding algorithms take a source image (src) and a threshold value (thresh) as input and produce an output image (dst) by comparing the pixel value at source pixel ( x , y ) to the threshold. If src ( x , y ) > thresh , then dst ( x , y ) is assigned some value. Otherwise dst ( x , y ) is assigned some other value.
In Binary Thresholding, in addition to the source image (src) and threshold value (thresh), it takes another input parameter called maximum value ( maxValue ). At each pixel location (x,y) it compares the pixel value src ( x , y ) to thresh. If src ( x , y ) is greater than thresh, it sets the value of the destination image pixel dst ( x , y ) to maxValue, otherwise it sets it to zero.

After thresholding I used the function cv.filter2D() to convolve a kernel with the image.The  kernel was defined using (kernel = numpy.ones((3, 3), numpy.float32) / 9) in our case. LPF (Low pass filter) helps in removing noise 
Filtering with the above kernel results in the following being performed: for each pixel, a 3x3 window is centered on this pixel, all pixels falling within this window are summed up, and the result is then divided by 9. This equates to computing the average of the pixel values inside that window. This operation is performed for all the pixels in the image to produce the output filtered image. 
After this the function cv2.medianBlur() takes a median of all the pixels under kernel area and central element is replaced with this median value. This is highly effective against salt-and-pepper noise in the images. In other filters, central element is a newly calculated value which may be a pixel value in the image or a new value. But in median blurring, central element is always replaced by some pixel value in the image. It reduces the noise effectively. 

Finally, the function cv2.dilate() which is a morphological operator is used. It is normally performed on binary images. It needs two inputs, one is our original image, the second one is called structuring element or kernel which decides the nature of operation. Here, a pixel element is ‘1’ if atleast one pixel under the kernel is ‘1’. So it increases the white region in the image or size of foreground object increases. Since noise is gone, they won’t come back, but our object area increases. It is also useful in joining broken parts of an object.
 
 
The user defined images to be analysed are firstly converted to grayscale ( gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY))
The color will vary from 0 to 255. The average color per row is calculated and then the average of these values is calculated to obtain a global average (using numpy.average).
We set the thresh value to be equal to the calculated average value - 1 and the maxValue to be equal to 255 in the above mentioned function.
After this step the function cv2.filter2D is used as described above.
The second and the third steps are used again after 4th sep is performed. But the thresh value is used is the average color obtained not the average color - 1.
After this median blur and dilate functions are used on the resulting image using the above mentioned kernel.



2.   How does the program detect the cells (as pixels)?
	It uses the function cv2.imread() to read the image and store as an ndarray of RGB, intensity values. You can access a pixel value by its row and column coordinates. For BGR image, it returns an array of Blue, Green, Red values. For grayscale image, just corresponding intensity is returned.
      
3.  How does the program fit the data (from the pixels ) and find the angle? 
    Linear Regression using the function LinearRegression().fit is performed on each thresholded (Binary image) frame of the cropped region containing the cell(Black pixels) to calculate the regression coefficients (slope and intercept). The slope is later used to find the angle the regression line makes with the positive x axis for each frame.
    The following steps are performed to find out the change in angle/frame and also assign a rotational direction (CCW or CW) (since slope to angle different for different quadrants can’t simply subtract. ) 

1. Difference in the consecutive values is taken to calculate the change in angle/frame  
    ith- (i-1)thlet us call this vector A.
2. Two more vectors are made by adding  180 to A i.e (Ai+ 180 ) and subtracting 180 i.e (Ai- 180) call them B , C respectively.
3. Another vector let’s say D is created whose ithelement is the Minimum of the absolute values of (Ai,Bi,Ci) . This will be our corrected value of change in angle.
4.  Vector E’s ith is  1 or -1 (IF( Ai/Di=1OR Bi/Di=1OR Ci/Di=1) then it is 1 otherwise -1 ), where -1 represents CW and +1 represents CCW.
5.  FPS can either be user defined or the program uses the function json.load and calculates the FPS value from the total elapsed time data from the metadata file present in the folder in which the stack of images of each frame of the rotating cell is present. It is used to calculate the change in angle/second which can further be used to calculate frequency using the formula 
      ω = 2πf.  (Program makes changes from degrees to radians wherever necessary).

 
4.   How does the program find the length of the cell, trace diameter, COM, axis of rotation (radius)? 
      
     Length of the cell is calculated by finding the diameter the circumcircle around the cell using the function cv2.minEnclosingCircle(). It is a circle which completely covers the object with minimum area. The diameter of the minimum enclosing circle is calculated for each frame and an average value gives us the length of the cell in cases where the standard deviation is high, mode is taken. 


   Trace Diameter is calculated by superimposing consecutive frames and using the function    cv2.findContours(). Contours can be explained simply as a curve joining all the continuous points (along the boundary), having the same color or intensity. The contours are a useful tool for shape analysis and object detection and recognition. We have thresholded binary images which in our case are inverted and then superimposed over one another we obtain a circular smudge of white pixels over a dark background. The findcontours function gives us the boundaries so that we can make use of the cv2.minEnclosingCircle() function again to find the diameter of the minimum enclosing circle.
      
     Centroid for each frame is found out by taking an average of the X coordinates, Y coordinates        of the pixels (Binary image) by using numpy. 
        xcom = numpy.mean(points[:, 0])
        ycom = numpy.mean(points[:, 1])
     The coordinates of the centroid are stored in a list as tuples. 

Radius, axis of rotation : 
We demonstrate 3 methods to find the center of Rotation of the system of points(pixels) using the above mentioned Centeriod data.
      We have our points in 2 list x and y.

a) Algebraic method
In this, we first find the mean of the points

x_m = mean(x), y_m = mean(y)

Next, we shift the origin of our coordinate system to x_m and y_m to get a reduced coordinate system
Calculation of the reduced coordinates
u = x - x_m
v = y - y_m 

To get the coordinates of CoM (uc, vc) in this reduced coordinate system, we solve the system of linear equation given below:

Suu * uc +  Suv * vc = (Suuu + Suvv)/2
Suv * uc +  Svv * vc = (Suuv + Svvv)/2
where,
Suv = sum(u * v)
Suu = sum(u**2)
Svv = sum(v**2)
Suuv = sum(u**2 * v)
Suvv = sum(u * v**2)
Suuu = sum(u**3)
Svvv = sum(v**3)

linalg.solve() function is used to solve these, giving (uc, vc) as the output.

Now to get back to our original coordinate system, we do the following.
xc_1 = x_m + uc
yc_1 = y_m + vc

Root mean square of the distance of all points from (xc_1, yc_1) gives us the radius R_1.

b) Least Square method
We use  optimize.leastsq() from the scipy library to get the coordinate of CoM.
It takes two arguments:
	1) Argument 1 -> A function to calculate the least square of the points the center
	2) An estimate of the center, calculated by taking the mean of x and y
	   center_estimate = x_m, y_m
	   where, x_m = mean(x)
	               y_m = mean(y)

The function optimize.leastsq(func, center_estimate) then outputs the desired coordinates (xc_2, yc_2)
Again, root mean square of the distance of all points from (xc_2, yc_2) gives us the radius R_2.

c) ODR method - Orthogonal Distance Regression
This is an advanced method to find the CoM which takes into account the errors present in the data itself. It is implemented in the scipy library.
First we have to instantiate ODR with your data, model and initial parameter estimates.
odr = odr.ODR(data, model, beta)
here, data is out data of points
      model is the model that we used to calculate the radius, i.e. least square value
      beta is the initial value [x_m, y_m, R_m]
      	where x_m = mean(x)
                             y_m = mean(y)
      	           R_m = initial guess of radius using (x_m, y_m)
We finally run the function to fit the data
output = odr.run()
which gives us the xc_3, yc_3, R_3 values.

We also plot the sum of square of differences to see which radius and CoM gives the best result.