description.txt

Jonathan's Photoelastic Disk Solver (TODO: Discuss a good name for the software) is a publicly available \cite{github} MATLAB implementation of the concept outlined in this manuscript. This program itself is split into three parts: Part 1, joDiskPrep.m, is a particle detector, neighbor finder and contact validator. Part 2, joDiskSolve.m, the actual photo elastic disk solver, and Part 3, joForceAdjMat.m, is a postprocessing and visualization tool.

Part 1, the preprocessor, takes an image or a series of images and detects circular particles in them, then identifies potential contacts by first coparing the surface distance between any two particles. If they are closer than the sum of their radii plus some tolerance theshold the photoelastic response in an area around the point of contact in both particles is evaluated and compared to a threshold value. If both particles show sufficient response near the contact point, the contact is declared valid. To get the best possible result it is necessary to adjust the contact threshold and the size of the evaluated area aorund the contact point for each specific set of images. A verbose option in the program allows the user to get a visualization of the detected contact network, plotted as an overlay to the camera image. When in doubt, the threshold should be set in a way that the preprocessor will detect more contacts than one think are valid. While this will result in increased computational load during the following fitting setps, the force solver is free to set the force of any contact to zero while it would not be allowed to create new contacts. Also several material- and experiental parameters are set here, such as the length scaling (particle size) and the photoelastic stress coefficient. The contact information and material parameters are then combined into a object oriented particle   based data structure. Finally, the combined set of data for all particles is saved into a single file for each image processed. 

Part 2, the solver, will process the files generated by the preprocessor and fits normal and tangential forces to all valid contacts, one particle at a time.
For each particle with a given number of contacts and their locations an initial guess for the magnitude of each contact force is calculated based on the gradiend squared method \cite{...}. This is done by first calculating total force on particle and then partiton that for all its contacts based on the sum of squares of the gradient in a small area aorund each contact point. An initial guess for the angles of attack $\alpha_i$ of each contact is generated by selecting them in such a way that with the given contact locations and magnitudes of forces force balance on the particle as a whole would be statisfied. This is and optional step that will dramaticaly improve the quality of the fit but might not be applicaple to all types of experiments. Then the fitting function is generated, that is a function that will generate a synthetic image of the photoelastic respone for a particle with the set contacts and free parameters for each contact force magnitude and angle of attack.  
This step is parallelized so the program can take advantage of mondern multi-core processors or be run on a high performance computer. Then a distance function $e$ is defined that gives the sum of squares of the difference between the experimental photoelastic response image $I_{ops}$ and its synthetic counterpart I_{synth}.

$$  e(f_{1..n}, \alpha_{1...n} ) = \sum_{x,y}{(I_{obs}(x,y)-I_{synth}(x,y,f_{1..n}, \alpha_{1...n}))^2} $$

This distance function is then fed into an optimizer (the default is the Levenberg-Marquard implementation contained in MATLAB's \emph{lsqnonlin}) and the $f_{1..n}, \alpha_{1...n}$ are varied to finde the minimum value for $d$. Note however that this will only find a local minimum in the $d(f,\alpha)$ landscape. Optionally, after each iteration of the optimizer the force balance condition for the particle is reapplied. As discussed above, this will greatly improve the accuracy of the fit as well as speed up the convergence. The fit results for $f$ and $\alpha$ and the residual value of $d$ are saved into the particle data structue. 

User defined parameters include detailed settings regarding the fitting algorithm and stopping criteria such as the maximum number of iterations of a tolerance value for the distance function.  Further included are and a scaling option to process a downsampled version of the photoelastic response image (to speed up processing) and an optional masking radius that can exclude the outer rim of a particle which can be prone to artefacts for example from residual stresses. In batch mode the program will process all valid preprocessing files found in the working directory. This allows it to run autonomously on a sperate computer, e.g. a high performance computing cluster.  
Finally, there is a verbose option available that allows to watch the fit for each particle converge graphically, that is, an the experimental image will be displayed next to an synthetic image based on the current iteration. 

Once all the particles of an input file are processed the results are saved into a single file again.

Part 3, the postprocessor reads the output files generated by the solver and then, for each contact, compares it to user defined quality criteria such as minimum and maximum values for the magnitude of the force or a maximum of the residual. Further, since each contact is actually fitted twice, one can evaluate is a contact obeys Newtwon's third law. If, for example due to image noise, the forces from each side of a contact point are detected slightly different one can then select to assign the average ob both forces to the contact or flag it for further inspection. If a contact passes all tests it's contact force is saved into an adjacency matrix representing the contact network in the current frame. Once all contacts have been processed the complete adjacency matrix is exported as well as a fully synthetic version of the force response image of the whole system. In verbose mode the fitted and validated contact network is drawn as an overlay on the experimental and synthetic force images.