Commit 1f346072 authored by Olivier Lantsoght's avatar Olivier Lantsoght
Browse files

[New] First cleaned code

# -*- coding: utf-8 -*-
Created on Mon Sep 4 11:20:14 2017
Example of script for force estimation based on photoelastic signal
@author: olantsoght
import PEGS as PEGS
# Filename and parameters
filenameMat ="Example/Polarized_end_cadre_solved.mat"
filenameFig ="Example/Polarized_end_cadre.JPG"
# Loading mat file and creating python equivalent structure
experiment = PEGS.load_matlab(
filenameMat, fitErrorCritera=0.42, g2Cal=0.5, brightfield=True,
verbose=True, pictureName=filenameFig, recomputeFit=True,
# Force identification
baseName = "Example/Polarized_end_cadre_solved"
print "First guess done"
experiment.forces_propagation(maxiter=5, verbose=True,
print "Propagation done"
print "Brute force done"
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
% Jonathan's Photoelastic Disk Solver
% Takes input from preprocessing script joDiskPrep.m as of 2016/09/27
% inspired from peDiskSolve by James Puckett (Phd-Thesis 2012)
% If you use this sovler please cite the follwoing paper
% K.E. Daniels, J. E. Kollmer & J. G. Puckett, "Photoelastic force measurements in granular materials", Rev. Sci. Inst. (201X)
% The generation of the synthetic force images is parallelized
% in a way that each row of the ouput image can be its own worker/thread
% it is usually limited by your number of processing cores though.
% Running this script may take a long time (up to one minute per particle
% on a 2016 desktop computer) so it makes sense to send this off to a high
% performance computer.
% Depending on your experimental data you can play around with the fit
% options, which might speed up processing considerabely.
% Last edit on 2016/09/28 by Jonathan Kollmer (
close all %housekeeping
clear all %housekeeping
%To do: correct G2Cal
G2Cal = 0.5;
%User defined values:
%which files are we processing ?
directory = '../Example/';
files = dir([directory, 'Polarized_end_cadre_preprocessing.mat']);
%how much of the particle diameter is used to fit the synthetic image
%(1 = use everything). Change this parameter only if the fit doesn't work
%correctly because of imaging artifacts at the particle boundary.
maskradius = 0.96;%
scaling = 1.0; %scale the image by this factor before doing the fit
brighfield = true;
%do we want to see each particle on screen while it is fitted ?
verbose = true;
%fit options: play around with these, depending on the quality of your
%images and the accuracy you want your result to have.
%Setting a good TolFun value can considerabely speed up your processing.
fitoptions = optimoptions('lsqnonlin','Algorithm','levenberg-marquardt','MaxIter',100,'MaxFunEvals',400,'TolFun',0.00001,'Display','final-detailed');
fitoptions_first = optimoptions('lsqnonlin','Algorithm','levenberg-marquardt','MaxIter',100,'MaxFunEvals',400,'TolFun',0.00001,'Display','final-detailed');
%fitoptions = optimoptions('lsqnonlin','Algorithm','trust-region-reflective','MaxIter',100,'MaxFunEvals',400,'TolFun',0.00001,'Display','final-detailed');
%There should be no need for user input below this line
nFrames = length(files); %how many files are we processing ?
for frame = 1:nFrames %loop over these frames
fileName = [directory,files(frame).name]; %which file/frame are we processing now ?
maskradius = maskradius / 2; %I did an unwise choice in naming this
load(fileName); %load the particle data file to process
pres=particle; %initialize a copy of the particle vector to mess around with
N = length(particle); %number of particles in this frame
display(['processing file ',fileName, 'containing ' ,num2str(N), 'particles']); %status indicator
for n=1:N
display(['fitting force(s) to particle ',num2str(n)]); %status indicator
if (particle(n).z > 0 )
z = particle(n).z;
fsigma = particle(n).fsigma;
rm = particle(n).rm;
g2 = particle(n).g2;
%This is the Camera Image we went to fit
template = particle(n).forceImage;
template = imresize(template,scaling);
%size of the force image
px = size(template,1);
%plot the experimental image
if verbose
%create a circular mask
template = histeq(template).*c_mask;
% To help the solver to get rid of false contact, we firstly %
% test all possible combination of active contact set %
% Number of possible combination
nbComb = zeros(1,z);
factZ = factorial(z);
for k=1:z
nbComb(k) = factZ/(factorial(k)*factorial(z-k));
% Creating combination
forcesIndex = zeros(sum(nbComb), z);
for k=1:z
forcesIndex(sum(nbComb(1:k))-nbComb(k)+1:sum(nbComb(1:k)),:) = [combnk(1:z,k), zeros(nbComb(k),z-k)];
% To store combination results
error = zeros(1,sum(nbComb));
error_immse = zeros(1,sum(nbComb));
error_ssim = zeros(1,sum(nbComb));
error_psnr = zeros(1,sum(nbComb));
param0= zeros(sum(nbComb), 2*z);
param = zeros(sum(nbComb), 2*z);
% Testing all combination
parfor comb = 1:sum(nbComb)
locZ = nnz(forcesIndex(comb,:));
% Uniform forces repartition, understimated
forces = (ones(1,locZ)*g2/G2Cal)/(locZ);
alphas = zeros(1,locZ);
beta = particle(n).betas(forcesIndex(comb,1:locZ));
p0 = [forces, alphas]
% Function computing intensity, "par" is [forces, alphas]
func = @(par) joForceImg(locZ, par(1:locZ),par(locZ+1:2*locZ), beta(1:locZ), fsigma, rm, px, 0);
% Error: difference between comptuted and experimental image
%err = @(par) psnr(func(par).*c_mask, template);
%err = @(par) immse(func(par).*c_mask, template);
%err = @(par) 1-ssim(func(par).*c_mask, template);
err = @(par) 1-ssim(histeq(func(par)).*c_mask, template);
%err = @(par) real(sum(sum( ( c_mask.*(template-func(par)).^2) )));
%BUG: for some reason I sometimes get imaginary results, this should not happen
% Solving
% Storing result
param0(comb,:) = [p0, zeros(1,2*z-2*locZ)];
param(comb,:) = [p, zeros(1,2*z-2*locZ)];
error(comb) = err(p);
error_immse(comb) = immse(func(p).*c_mask, template);
error_ssim(comb) = ssim(func(p).*c_mask, template)
error_psnr(comb) = psnr(func(p).*c_mask, template);
display(['Particle ',num2str(n),': combination ', num2str(comb), ' over ', num2str(sum(nbComb)), ' done']); %status indicator
% Retrieve best candidate and run again optim
[~,bestIndex] = min(error);
locZ = nnz(forcesIndex(bestIndex,:));
beta = particle(n).betas(forcesIndex(bestIndex,1:locZ));
p0 = param(bestIndex, 1:2*locZ);
% Function computing intensity, "par" is [forces, alphas]
func = @(par) joForceImg(locZ, par(1:locZ),par(locZ+1:2*locZ), beta(1:locZ), fsigma, rm, px, verbose);
% Error: difference between comptuted and experimental image
%err = @(par) psnr(func(par).*c_mask, template);
%err = @(par) immse(func(par).*c_mask, template);
%err = @(par) 1-ssim(func(par).*c_mask, template);
err = @(par) 1-ssim(histeq(func(par)).*c_mask, template);
%err = @(par) real(sum(sum( ( c_mask.*(template-func(par)).^2) )));
% Solving
%get back the result from fitting
forces = p(1:locZ);
alphas = p(locZ+1:2*locZ);
fitError = err(p);
%generate an image with the fitted parameters
imgFit = joForceImg(locZ, forces, alphas, beta, fsigma, rm, px*(1/scaling), verbose);
%store the new information into a new particle vector
pres(n).z = locZ;
pres(n).f = g2/G2Cal;
pres(n).forces = forces;
pres(n).betas = beta;
pres(n).alphas = alphas;
pres(n).neighbours = particle(n).neighbours(forcesIndex(bestIndex,1:locZ));
% TO DO : removing antagonist contact? Not a god idea if it did
% not well converged
pres(n).contactG2s = particle(n).contactG2s(forcesIndex(bestIndex,1:locZ));
%pres(n).contactIs = particle(n).contactIs(forcesIndex(bestIndex,1:locZ));
% Not always exist
pres(n).fitError = fitError;
pres(n).synthImg = imgFit;
%save the result
time = toc
%Computes the optical response of any given point in a loaded photoelastic disk
%Originally from James G. Puckett and Karen E. Daniels. "Equilibrating temperaturelike variables in jammed granular subsystems." Physical Review Letters 110 058001 (2013)
%Ported to Matlab by Jonathan E. Kollmer, 2016-05-06
%result = photoelastic response (intensity)
%xxi = x-position at which the stress is computed (in meters)
%xxy = y-position at which the stress is computed (in meters)
%z = coordination number
%f list of forces (magintude) acting on the particle
%alpha list of angles of the forces, see diagram in J.Puckett, PhD-Thesis, North Carolina State University
%beta list of angles of the forces, see diagram in J.Puckett, PhD-Thesis, North Carolina State University
%fsigma photoelastic stress coefficient
%rm radius of the disk (in meters)
function result = StressEngine(xxi, xxj, z, f, alpha, beta, fsigma, rm)
beta = -beta+pi/2;
pioverfsigma = pi/(fsigma);
oneoverpirm = 1./(pi*rm);
twooverpi = 2./pi;
%double aa,a,b,b2,x1,y1,x2,y2,ch0,ch1,chn,costh1,r10,r11,r1n,th1,s1, s2,sr,th; %temporary
for k=1:z
b = beta(k)+pi/2; %rotate pi/2, to match input image
a = alpha(k);
if (a<0)
b2 = b+(pi+2*a);
b2 = b-(pi-2*a);
x1 = rm*sin(b);
y1 = rm*cos(b);
x2 = rm*sin(b2);
y2 = rm*cos(b2);
ch0 = x2-x1; % chord x
ch1 = y2-y1; % chord y
chn = sqrt(ch0^2+ch1^2);
ch0 = ch0/chn; % normalize chord x
ch1 = ch1/chn; % normalize chord y
r10 = xxi - x1; % r vector x coord
r11 =-xxj - y1; % r vector y coord
r1n = sqrt(r10^2+r11^2);
costh1 = (r10*ch0+r11*ch1)/r1n;
if (r11*ch0>r10*ch1) % important!
signth = 1;
signth = -1;
th1 = signth*acos(costh1); %faster than cos(asin(stuff));
s2 = -(f(k)*oneoverpirm)*(-sin(a)); %uniform compression from boundary, sin(pi-a) = -sin(a)
s1 = -(f(k)*twooverpi)/(r1n)*costh1;
sr = s1-s2;
th = th1-beta(k)-alpha(k); %rotate coordinates
sigmaxx = sigmaxx + sr*((sin(th))^2);
sigmayy = sigmayy + sr*((cos(th))^2);
sigmaxy = sigmaxy + 0.5*s1*(sin(2*th));
end %end of k loop over beta, alpha
aa = real(sqrt((sigmaxx-sigmayy)^2+4*(sigmaxy)^2));
result = (cos(pioverfsigma*aa))^2; %wrap = sin(2*pi*t/fsigma*a).^2;
if (isnan(result)) %for some reason result sometimes gets to be NAN, not sure why this happens
result = 0; %temproary fix is to set it zero then
if (result<0) %for some reason, not sure why this happens
result = 0; %temproary fix is to set it zero then
function [z, r, residual] = fitcircle(x, varargin)
%FITCIRCLE least squares circle fit
% [Z, R] = FITCIRCLE(X) fits a circle to the N points in X minimising
% geometric error (sum of squared distances from the points to the fitted
% circle) using nonlinear least squares (Gauss Newton)
% Input
% X : 2xN array of N 2D points, with N >= 3
% Output
% Z : center of the fitted circle
% R : radius of the fitted circle
% [Z, R] = FITCIRCLE(X, 'linear') fits a circle using linear least
% squares minimising the algebraic error (residual from fitting system
% of the form ax'x + b'x + c = 0)
% [Z, R] = FITCIRCLE(X, Property, Value, ...) allows parameters to be
% passed to the internal Gauss Newton method. Property names can be
% supplied as any unambiguous contraction of the property name and are
% case insensitive, e.g. FITCIRCLE(X, 't', 1e-4) is equivalent to
% FITCIRCLE(X, 'tol', 1e-4). Valid properties are:
% Property: Value:
% --------------------------------
% maxits positive integer, default 100
% Sets the maximum number of iterations of the Gauss Newton
% method
% tol positive constant, default 1e-5
% Gauss Newton converges when the relative change in the solution
% is less than tol
% [X, R, RES] = fitcircle(...) returns the 2 norm of the residual from
% the least squares fit
% Example:
% x = [1 2 5 7 9 3; 7 6 8 7 5 7];
% % Get linear least squares fit
% [zl, rl] = fitcircle(x, 'linear')
% % Get true best fit
% [z, r] = fitcircle(x)
% Reference: "Least-squares fitting of circles and ellipses", W. Gander,
% G. Golub, R. Strebel - BIT Numerical Mathematics, 1994, Springer
% This implementation copyright Richard Brown, 2007, but is freely
% available to copy, use, or modify as long as this line is maintained
error(nargchk(1, 5, nargin, 'struct'))
% Default parameters for Gauss Newton minimisation
params.maxits = 100;
params.tol = 1e-5;
% Check x and get user supplied parameters
[x, fNonlinear, params] = parseinputs(x, params, varargin{:});
% Convenience variables
m = size(x, 2);
x1 = x(1, :)';
x2 = x(2, :)';
% 1) Compute best fit w.r.t. algebraic error using linear least squares
% Circle is represented as a matrix quadratic form
% ax'x + b'x + c = 0
% Linear least squares estimate found by minimising Bu = 0 s.t. norm(u) = 1
% where u = [a; b; c]
% Form the coefficient matrix
B = [x1.^2 + x2.^2, x1, x2, ones(m, 1)];
% Least squares estimate is right singular vector corresp. to smallest
% singular value of B
[U, S, V] = svd(B);
u = V(:, 4);
% For clarity, set the quadratic form variables
a = u(1);
b = u(2:3);
c = u(4);
% Convert to centre/radius
z = -b / (2*a);
r = sqrt((norm(b)/(2*a))^2 - c/a);
% 2) Nonlinear refinement to miminise geometric error, and compute residual
if fNonlinear
[z, r, residual] = fitcircle_geometric(x, z, r);
residual = norm(B * u);
function [z, r, residual] = fitcircle_geometric(x, z0, r0)
% Use a simple Gauss Newton method to minimize the geometric error
fConverged = false;
% Set initial u
u = [z0; r0];
% Delta is the norm of current step, scaled by the norm of u
delta = inf;
nIts = 0;
for nIts = 1:params.maxits
% Find the function and Jacobian
[f, J] = sys(u);
% Solve for the step and update u
h = -J \ f;
u = u + h;
% Check for convergence
delta = norm(h, inf) / norm(u, inf);
if delta < params.tol
fConverged = true;
if ~fConverged
warning('fitcircle:FailureToConverge', ...
'Gauss Newton iteration failed to converge');
z = u(1:2);
r = u(3);
f = sys(u);
residual = norm(f);
function [f, J] = sys(u)
%SYS Nonlinear system to be minimised - the objective
%function is the distance to each point from the fitted circle
%contained in u
% Objective function
f = (sqrt(sum((repmat(u(1:2), 1, m) - x).^2)) - u(3))';
% Jacobian
denom = sqrt( (u(1) - x1).^2 + (u(2) - x2).^2 );
J = [(u(1) - x1) ./ denom, (u(2) - x2) ./ denom, repmat(-1, m, 1)];
end % sys
end % fitcircle_geometric
end % fitcircle
function [x, fNonlinear, params] = parseinputs(x, params, varargin)
% Make sure x is 2xN where N > 3
if size(x, 2) == 2
x = x';
if size(x, 1) ~= 2
error('fitcircle:InvalidDimension', ...
'Input matrix must be two dimensional')
if size(x, 2) < 3
error('fitcircle:InsufficientPoints', ...
'At least 3 points required to compute fit')
% determine whether we are measuring geometric error (nonlinear), or
% algebraic error (linear)
fNonlinear = true;
switch length(varargin)
% No arguments means a nonlinear least squares with defaul parameters
case 0
% One argument can only be 'linear', specifying linear least squares
case 1
if strncmpi(varargin{1}, 'linear', length(varargin{1}))
fNonlinear = false;
error('fitcircle:UnknownOption', 'Unknown Option')
% Otherwise we're left with user supplied parameters for Gauss Newton
if rem(length(varargin), 2) ~= 0
error('fitcircle:propertyValueNotPair', ...
'Additional arguments must take the form of Property/Value pairs');
% Cell array of valid property names
properties = {'maxits', 'tol'};
while length(varargin) ~= 0
property = varargin{1};
value = varargin{2};
% If the property has been supplied in a shortened form, lengthen it
iProperty = find(strncmpi(property, properties, length(property)));
if isempty(iProperty)
error('fitcircle:UnkownProperty', 'Unknown Property');
elseif length(iProperty) > 1
error('fitcircle:AmbiguousProperty', ...
'Supplied shortened property name is ambiguous');
% Expand property to its full name
property = properties{iProperty};
switch property
case 'maxits'
if value <= 0
error('fitcircle:InvalidMaxits', ...
'maxits must be an integer greater than 0')
params.maxits = value;
case 'tol'
if value <= 0
error('fitcircle:InvalidTol', ...
'tol must be a positive real number')
params.tol = value;
end % switch property
varargin(1:2) = [];
end % while
end % switch length(varargin)
end %parseinputs
\ No newline at end of file
%This function creates a synthetic photoelastic response image
%last change on 2016/09/27 by Jonathan Kollmer
% %test values
% z=2; %Number of contacts this particle has
% f = [0.6 0.6]; %Absolute forces on this particle
% alpha = [0 0]; %Alpha contact angles on this particle
% beta = [0 -pi]; %Beta contact angles on this particle
% fsigma = 100; %Photoelastic stress coefficient of this particle
% rm = 0.00816516; %Particle radius in meters
% px = 100; %Return image size in pixels
function img = joForceImg (z, f, alpha, beta, fsigma, rm, px, verbose)
%make sure the forces are balanced
%[alpha,f] = forceBalance2(f,alpha,beta);
%create an empty placeholder image for our result/return value
img = zeros(px);
%Create a scale that maps the diameter of the particle onto the image size
xx = linspace(-rm, rm, px);
parfor (x=1:px) %loop throgh image width
xRow=ones(px,1); %placeholder for the current row beeing prcessed, parallel for makes it necessary to split up the result data and later combine it
for y=1:px %loop throgh image height
if ((xx(x)^2+xx(y)^2)<=rm^2) %check if we are actually inside
xRow(y) = StressEngine(xx(x), xx(y), z, f, alpha, beta, fsigma, rm); %call the StressEngine to compute the photoelastic response at each pixel
img(x,:)=xRow; %consolidate processed row into the output image, necessary data shuffling to use parallel for