This adds manually designed features to the OpenSeeData.

This commit is contained in:
Emiliana 2020-03-01 14:23:56 +01:00
parent 7706372a7a
commit 04d7f5b8b3
7 changed files with 421 additions and 49 deletions

19
Licenses/Remedian.txt Normal file
View file

@ -0,0 +1,19 @@
Copyright (c) 2016, Tim Menzies tim@menzies.us, MIT license v2.
Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
the Software without restriction, including without limitation the rights to
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
of the Software, and to permit persons to whom the Software is furnished to do
so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

View file

@ -16,7 +16,7 @@ public class OpenSee : MonoBehaviour {
public int listenPort = 11573; public int listenPort = 11573;
private const int nPoints = 68; private const int nPoints = 68;
private const int packetFrameSize = 8 + 4 + 2 * 4 + 2 * 4 + 1 + 4 + 3 * 4 + 3 * 4 + 4 * 4 + 4 * 68 + 4 * 2 * 68 + 4 * 3 * 70; private const int packetFrameSize = 8 + 4 + 2 * 4 + 2 * 4 + 1 + 4 + 3 * 4 + 3 * 4 + 4 * 4 + 4 * 68 + 4 * 2 * 68 + 4 * 3 * 70 + 4 * 14;
[Header("Tracking data")] [Header("Tracking data")]
[Tooltip("This is an informational property that tells you how many packets have been received")] [Tooltip("This is an informational property that tells you how many packets have been received")]
@ -58,6 +58,40 @@ public class OpenSee : MonoBehaviour {
public Vector2[] points; public Vector2[] points;
[Tooltip("These are 3D points estimated from the 2D points. The should be rotation and translation compensated. There are 70 points with guesses for the eyeball center positions being added at the end of the 68 2D points.")] [Tooltip("These are 3D points estimated from the 2D points. The should be rotation and translation compensated. There are 70 points with guesses for the eyeball center positions being added at the end of the 68 2D points.")]
public Vector3[] points3D; public Vector3[] points3D;
[Tooltip("This field contains a number of action unit like features.")]
public OpenSeeFeatures features;
[System.Serializable]
public class OpenSeeFeatures {
[Tooltip("This field indicates whether the left eye is opened(0) or closed (-1). A value of 1 means open wider than normal.")]
public float EyeLeft;
[Tooltip("This field indicates whether the right eye is opened(0) or closed (-1). A value of 1 means open wider than normal.")]
public float EyeRight;
[Tooltip("This field indicates how steep the left eyebrow is, compared to the median steepness.")]
public float EyebrowSteepnessLeft;
[Tooltip("This field indicates how far up or down the left eyebrow is, compared to its median position.")]
public float EyebrowUpDownLeft;
[Tooltip("This field indicates how quirked the left eyebrow is, compared to its median quirk.")]
public float EyebrowQuirkLeft;
[Tooltip("This field indicates how steep the right eyebrow is, compared to the average steepness.")]
public float EyebrowSteepnessRight;
[Tooltip("This field indicates how far up or down the right eyebrow is, compared to its median position.")]
public float EyebrowUpDownRight;
[Tooltip("This field indicates how quirked the right eyebrow is, compared to its median quirk.")]
public float EyebrowQuirkRight;
[Tooltip("This field indicates how far up or down the left mouth corner is, compared to its median position.")]
public float MouthCornerUpDownLeft;
[Tooltip("This field indicates how far in or out the left mouth corner is, compared to its median position.")]
public float MouthCornerInOutLeft;
[Tooltip("This field indicates how far up or down the right mouth corner is, compared to its median position.")]
public float MouthCornerUpDownRight;
[Tooltip("This field indicates how far in or out the right mouth corner is, compared to its median position.")]
public float MouthCornerInOutRight;
[Tooltip("This field indicates how open or closed the mouth is, compared to its median pose.")]
public float MouthOpen;
[Tooltip("This field indicates how wide the mouth is, compared to its median pose.")]
public float MouthWide;
}
public OpenSeeData() { public OpenSeeData() {
confidence = new float[nPoints]; confidence = new float[nPoints];
@ -139,6 +173,22 @@ public class OpenSee : MonoBehaviour {
rightGaze = Quaternion.LookRotation(swapX(points3D[66]) - swapX(points3D[68])) * Quaternion.AngleAxis(180, Vector3.right) * Quaternion.AngleAxis(180, Vector3.forward); rightGaze = Quaternion.LookRotation(swapX(points3D[66]) - swapX(points3D[68])) * Quaternion.AngleAxis(180, Vector3.right) * Quaternion.AngleAxis(180, Vector3.forward);
leftGaze = Quaternion.LookRotation(swapX(points3D[67]) - swapX(points3D[69])) * Quaternion.AngleAxis(180, Vector3.right) * Quaternion.AngleAxis(180, Vector3.forward); leftGaze = Quaternion.LookRotation(swapX(points3D[67]) - swapX(points3D[69])) * Quaternion.AngleAxis(180, Vector3.right) * Quaternion.AngleAxis(180, Vector3.forward);
features = new OpenSeeFeatures();
features.EyeLeft = readFloat(b, ref o);
features.EyeRight = readFloat(b, ref o);
features.EyebrowSteepnessLeft = readFloat(b, ref o);
features.EyebrowUpDownLeft = readFloat(b, ref o);
features.EyebrowQuirkLeft = readFloat(b, ref o);
features.EyebrowSteepnessRight = readFloat(b, ref o);
features.EyebrowUpDownRight = readFloat(b, ref o);
features.EyebrowQuirkRight = readFloat(b, ref o);
features.MouthCornerUpDownLeft = readFloat(b, ref o);
features.MouthCornerInOutLeft = readFloat(b, ref o);
features.MouthCornerUpDownRight = readFloat(b, ref o);
features.MouthCornerInOutRight = readFloat(b, ref o);
features.MouthOpen = readFloat(b, ref o);
features.MouthWide = readFloat(b, ref o);
} }
} }

View file

@ -27,6 +27,8 @@ public class OpenSeeExpression : MonoBehaviour
public string filename = ""; public string filename = "";
[Tooltip("When enabled, training will be done in such a way that the model can output probabilities. This is a lot slower!")] [Tooltip("When enabled, training will be done in such a way that the model can output probabilities. This is a lot slower!")]
public bool enableProbabilityTraining = false; public bool enableProbabilityTraining = false;
[Tooltip("This is the SVM training C value.")]
public float C = 2f;
[Header("Toggles")] [Header("Toggles")]
[Tooltip("When enabled, calibration data will be collected from the given OpenSee component.")] [Tooltip("When enabled, calibration data will be collected from the given OpenSee component.")]
public bool recording = false; public bool recording = false;
@ -88,6 +90,8 @@ public class OpenSeeExpression : MonoBehaviour
public bool pointsLipUpper = true; public bool pointsLipUpper = true;
[Tooltip("When enabled, the points of the lower lip will be used to determine the expression.")] [Tooltip("When enabled, the points of the lower lip will be used to determine the expression.")]
public bool pointsLipLower = true; public bool pointsLipLower = true;
[Tooltip("When enabled, certain manually designed features will be used to determine the expression.")]
public bool features = true;
} }
[Serializable] [Serializable]
private class OpenSeeExpressionRepresentation { private class OpenSeeExpressionRepresentation {
@ -96,6 +100,7 @@ public class OpenSeeExpression : MonoBehaviour
private string[] classLabels = null; private string[] classLabels = null;
private int[] indices = null; private int[] indices = null;
private PointSelection pointSelection; private PointSelection pointSelection;
private bool newModel = false;
static public void LoadSerialized(byte[] modelBytes, out Dictionary<string, List<float[]>> expressions, out SVMModel model, out string[] classLabels, out int[] indices, out PointSelection pointSelection) { static public void LoadSerialized(byte[] modelBytes, out Dictionary<string, List<float[]>> expressions, out SVMModel model, out string[] classLabels, out int[] indices, out PointSelection pointSelection) {
IFormatter formatter = new BinaryFormatter(); IFormatter formatter = new BinaryFormatter();
@ -117,6 +122,8 @@ public class OpenSeeExpression : MonoBehaviour
} }
if (pointSelection == null) if (pointSelection == null)
pointSelection = new PointSelection(); pointSelection = new PointSelection();
if (!oser.newModel)
pointSelection.features = true;
} }
static public byte[] ToSerialized(Dictionary<string, List<float[]>> expressions, SVMModel model, string[] classLabels, int[] indices, PointSelection pointSelection) { static public byte[] ToSerialized(Dictionary<string, List<float[]>> expressions, SVMModel model, string[] classLabels, int[] indices, PointSelection pointSelection) {
@ -126,6 +133,7 @@ public class OpenSeeExpression : MonoBehaviour
oser.classLabels = classLabels; oser.classLabels = classLabels;
oser.indices = indices; oser.indices = indices;
oser.pointSelection = pointSelection; oser.pointSelection = pointSelection;
oser.newModel = true;
IFormatter formatter = new BinaryFormatter(); IFormatter formatter = new BinaryFormatter();
MemoryStream memoryStream = new MemoryStream(); MemoryStream memoryStream = new MemoryStream();
@ -141,9 +149,9 @@ public class OpenSeeExpression : MonoBehaviour
private Dictionary<string, List<float[]>> expressions; private Dictionary<string, List<float[]>> expressions;
private SVMModel model = null; private SVMModel model = null;
private string[] classLabels = null; private string[] classLabels = null;
private int maxSamples = 1000; // Allows 10 expressions without exceeding 10000 rows in total for SVM training, which is a commonly seen upper limit. private int maxSamples = 1666;
// rightEyeOpen, leftEyeOpen, translation, rawQuaternion, rawEuler, confidence, points/(width, height), points3D // rightEyeOpen, leftEyeOpen, translation, rawQuaternion, rawEuler, confidence, points/(width, height), points3D, features
private int colsFull = 1 + 1 + 3 + 4 + 3 + /*66 + 2 * 66 +*/ 3 * 66; private int colsFull = 1 + 1 + 3 + 4 + 3 + /*66 + 2 * 66 +*/ 3 * 66 + 14;
private int colsBase = 1 + 1 + 3 + 4 + 3; private int colsBase = 1 + 1 + 3 + 4 + 3;
private int cols; private int cols;
private double lastCapture = 0.0; private double lastCapture = 0.0;
@ -235,6 +243,11 @@ public class OpenSeeExpression : MonoBehaviour
indexList.Add(colsBase + i * 3 + 2); indexList.Add(colsBase + i * 3 + 2);
cols += 3; cols += 3;
} }
if (pointSelection.features)
for (int i = 0; i < 14; i++) {
indexList.Add(colsBase + 66 * 3 + i);
cols++;
}
indices = indexList.ToArray(); indices = indexList.ToArray();
} }
@ -294,10 +307,26 @@ public class OpenSeeExpression : MonoBehaviour
data[12 + /*66 + 2 * 66 +*/ 3 * i + 1] = t.points3D[i].y; data[12 + /*66 + 2 * 66 +*/ 3 * i + 1] = t.points3D[i].y;
data[12 + /*66 + 2 * 66 +*/ 3 * i + 2] = t.points3D[i].z; data[12 + /*66 + 2 * 66 +*/ 3 * i + 2] = t.points3D[i].z;
} }
data[12 + 66 * 3 + 0] = t.features.EyeLeft;
data[12 + 66 * 3 + 1] = t.features.EyeRight;
data[12 + 66 * 3 + 2] = t.features.EyebrowSteepnessLeft;
data[12 + 66 * 3 + 3] = t.features.EyebrowUpDownLeft;
data[12 + 66 * 3 + 4] = t.features.EyebrowQuirkLeft;
data[12 + 66 * 3 + 5] = t.features.EyebrowSteepnessRight;
data[12 + 66 * 3 + 6] = t.features.EyebrowUpDownRight;
data[12 + 66 * 3 + 7] = t.features.EyebrowQuirkRight;
data[12 + 66 * 3 + 8] = t.features.MouthCornerUpDownLeft;
data[12 + 66 * 3 + 9] = t.features.MouthCornerInOutLeft;
data[12 + 66 * 3 + 10] = t.features.MouthCornerUpDownRight;
data[12 + 66 * 3 + 11] = t.features.MouthCornerInOutRight;
data[12 + 66 * 3 + 12] = t.features.MouthOpen;
data[12 + 66 * 3 + 13] = t.features.MouthWide;
return data; return data;
} }
public bool TrainModel() { public bool TrainModel() {
Debug.Log("Training model");
Debug.Log("---------------------");
train = false; train = false;
if (openSee == null) { if (openSee == null) {
ResetInfo(); ResetInfo();
@ -310,7 +339,8 @@ public class OpenSeeExpression : MonoBehaviour
List<string> accuracyWarnings = new List<string>(); List<string> accuracyWarnings = new List<string>();
int samples = 0; int samples = 0;
foreach (string key in expressions.Keys) { foreach (string key in expressions.Keys) {
List<float[]> list = expressions[key]; List<float[]> list = new List<float[]>(expressions[key]);
list.RemoveAll(x => x.Length != colsFull);
if (list != null && list.Count == maxSamples) { if (list != null && list.Count == maxSamples) {
keys.Add(key); keys.Add(key);
samples += maxSamples; samples += maxSamples;
@ -337,7 +367,7 @@ public class OpenSeeExpression : MonoBehaviour
classLabels = keys.ToArray(); classLabels = keys.ToArray();
int train_split = maxSamples * 3 / 4; int train_split = maxSamples * 3 / 4;
int test_split = maxSamples - train_split; int test_split = maxSamples - train_split;
int rows_train = classes * train_split; int rows_train = (classes + 2) * train_split;
int rows_test = classes * test_split; int rows_test = classes * test_split;
float[] X_train = new float[rows_train * cols]; float[] X_train = new float[rows_train * cols];
float[] X_test = new float[rows_test * cols]; float[] X_test = new float[rows_test * cols];
@ -347,7 +377,8 @@ public class OpenSeeExpression : MonoBehaviour
int i_test = 0; int i_test = 0;
System.Random rnd = new System.Random(); System.Random rnd = new System.Random();
for (int i = 0; i < classes; i++) { for (int i = 0; i < classes; i++) {
List<float[]> list = expressions[keys[i]]; List<float[]> list = new List<float[]>(expressions[keys[i]]);
list.RemoveAll(x => x.Length != colsFull);
int local_train_split = list.Count * 3 / 4; int local_train_split = list.Count * 3 / 4;
test_split = list.Count - local_train_split; test_split = list.Count - local_train_split;
for (int j = list.Count; j > 1;) { for (int j = list.Count; j > 1;) {
@ -357,16 +388,28 @@ public class OpenSeeExpression : MonoBehaviour
list[k] = list[j]; list[k] = list[j];
list[j]= tmp; list[j]= tmp;
} }
for (int j = 0; j < train_split; j++) { int train_split_current = train_split;
if (classLabels[i] == "neutral")
train_split_current *= 3;
for (int j = 0; j < train_split_current; j++) {
float factor = 1f;
float adder = 0f;
if (j > train_split || j > list.Count) {
factor = UnityEngine.Random.Range(0.98f, 1.02f);
adder = UnityEngine.Random.Range(-0.02f, 0.02f);
}
for (int k = 0; k < cols; k++) { for (int k = 0; k < cols; k++) {
X_train[i_train * cols + k] = list[j % local_train_split][indices[k]]; float v = list[j % local_train_split][indices[k]] * factor + adder;
X_train[i_train * cols + k] = v;
} }
y_train[i_train] = i; y_train[i_train] = i;
i_train++; i_train++;
} }
for (int j = local_train_split; j < local_train_split + test_split; j++) { for (int j = local_train_split; j < local_train_split + test_split; j++) {
for (int k = 0; k < cols; k++) for (int k = 0; k < cols; k++) {
X_test[i_test * cols + k] = list[j][indices[k]]; float v = list[j][indices[k]];
X_test[i_test * cols + k] = v;
}
y_test[i_test] = i; y_test[i_test] = i;
i_test++; i_test++;
} }
@ -374,7 +417,7 @@ public class OpenSeeExpression : MonoBehaviour
int probability = 0; int probability = 0;
if (enableProbabilityTraining) if (enableProbabilityTraining)
probability = 1; probability = 1;
model.TrainModel(X_train, y_train, rows_train, cols, probability); model.TrainModel(X_train, y_train, rows_train, cols, probability, C);
confusionMatrix = model.ConfusionMatrix(X_test, y_test, i_test, out accuracy); confusionMatrix = model.ConfusionMatrix(X_test, y_test, i_test, out accuracy);
confusionMatrixString = SVMModel.FormatMatrix(confusionMatrix, classLabels); confusionMatrixString = SVMModel.FormatMatrix(confusionMatrix, classLabels);
for (int label = 0; label < classes; label++) { for (int label = 0; label < classes; label++) {

View file

@ -118,10 +118,10 @@ public class SVMModel {
} }
// features is a rows*cols long float array with rows written one after another. labels is a rows long float array with the corresponding classes as integral numbers from 0 to num_classes-1, where all classes have to exist. Training on more than 10000 rows is disabled. // features is a rows*cols long float array with rows written one after another. labels is a rows long float array with the corresponding classes as integral numbers from 0 to num_classes-1, where all classes have to exist. Training on more than 10000 rows is disabled.
public bool TrainModel(float[] features, float[] labels, int rows, int cols, int probability) { public bool TrainModel(float[] features, float[] labels, int rows, int cols, int probability, float C) {
DestroyModel(); DestroyModel();
if (rows > 10000) //if (rows > 10000)
return false; // return false;
if (cols < 1 || rows < 1 || features.Length < rows * cols || labels.Length < rows) if (cols < 1 || rows < 1 || features.Length < rows * cols || labels.Length < rows)
return false; return false;
@ -143,7 +143,7 @@ public class SVMModel {
maxClasses = max + 1; maxClasses = max + 1;
this.cols = cols; this.cols = cols;
model = trainModel(features, labels, rows, cols, maxClasses, probability, 2f); model = trainModel(features, labels, rows, cols, maxClasses, probability, C);
haveModel = true; haveModel = true;
return true; return true;

View file

@ -79,6 +79,8 @@ tracking_time = 0.0
tracking_frames = 0 tracking_frames = 0
frame_count = 0 frame_count = 0
features = ["eye_l", "eye_r", "eyebrow_steepness_l", "eyebrow_updown_l", "eyebrow_quirk_l", "eyebrow_steepness_r", "eyebrow_updown_r", "eyebrow_quirk_r", "mouth_corner_updown_l", "mouth_corner_inout_l", "mouth_corner_updown_r", "mouth_corner_inout_r", "mouth_open", "mouth_wide"]
if args.log_data != "": if args.log_data != "":
log = open(args.log_data, "w") log = open(args.log_data, "w")
log.write("Frame,Time,Width,Height,FPS,Face,FaceID,RightOpen,LeftOpen,AverageConfidence,Success3D,PnPError,RotationQuat.X,RotationQuat.Y,RotationQuat.Z,RotationQuat.W,Euler.X,Euler.Y,Euler.Z,RVec.X,RVec.Y,RVec.Z,TVec.X,TVec.Y,TVec.Z") log.write("Frame,Time,Width,Height,FPS,Face,FaceID,RightOpen,LeftOpen,AverageConfidence,Success3D,PnPError,RotationQuat.X,RotationQuat.Y,RotationQuat.Z,RotationQuat.W,Euler.X,Euler.Y,Euler.Z,RVec.X,RVec.Y,RVec.Z,TVec.X,TVec.Y,TVec.Z")
@ -86,6 +88,8 @@ if args.log_data != "":
log.write(f",Landmark[{i}].X,Landmark[{i}].Y,Landmark[{i}].Confidence") log.write(f",Landmark[{i}].X,Landmark[{i}].Y,Landmark[{i}].Confidence")
for i in range(66): for i in range(66):
log.write(f",Point3D[{i}].X,Point3D[{i}].Y,Point3D[{i}].Z") log.write(f",Point3D[{i}].X,Point3D[{i}].Y,Point3D[{i}].Z")
for feature in features:
log.write(f",{feature}")
log.write("\r\n") log.write("\r\n")
log.flush() log.flush()
@ -120,8 +124,10 @@ try:
for face_num, f in enumerate(faces): for face_num, f in enumerate(faces):
f = copy.copy(f) f = copy.copy(f)
f.id += args.face_id_offset f.id += args.face_id_offset
right_state = "O" if f.eye_blink[0] > 0.8 else "-" if f.eye_blink is None:
left_state = "O" if f.eye_blink[1] > 0.8 else "-" f.eye_blink = [1, 1]
right_state = "O" if f.eye_blink[0] > 0.30 else "-"
left_state = "O" if f.eye_blink[1] > 0.30 else "-"
if args.silent == 0: if args.silent == 0:
print(f"Confidence[{f.id}]: {f.conf:.4f} / 3D fitting error: {f.pnp_error:.4f} / Eyes: {left_state}, {right_state}") print(f"Confidence[{f.id}]: {f.conf:.4f} / 3D fitting error: {f.pnp_error:.4f} / Eyes: {left_state}, {right_state}")
detected = True detected = True
@ -154,9 +160,9 @@ try:
packet.extend(bytearray(struct.pack("f", x))) packet.extend(bytearray(struct.pack("f", x)))
if not log is None: if not log is None:
log.write(f",{y},{x},{c}") log.write(f",{y},{x},{c}")
if pt_num == 66 and f.eye_blink[0] < 0.5: if pt_num == 66 and f.eye_blink[0] < 0.30:
continue continue
if pt_num == 67 and f.eye_blink[0] < 0.5: if pt_num == 67 and f.eye_blink[1] < 0.30:
continue continue
x = int(x + 0.5) x = int(x + 0.5)
y = int(y + 0.5) y = int(y + 0.5)
@ -179,9 +185,9 @@ try:
x -= 1 x -= 1
if not (x < 0 or y < 0 or x >= height or y >= width): if not (x < 0 or y < 0 or x >= height or y >= width):
frame[int(x), int(y)] = color frame[int(x), int(y)] = color
if args.pnp_points != 0 and (args.visualize != 0 or not out is None): if args.pnp_points != 0 and (args.visualize != 0 or not out is None) and f.rotation is not None:
if args.pnp_points > 1: if args.pnp_points > 1:
projected = cv2.projectPoints(f.face_3d, f.rotation, f.translation, tracker.camera, tracker.dist_coeffs) projected = cv2.projectPoints(f.face_3d[0:66], f.rotation, f.translation, tracker.camera, tracker.dist_coeffs)
else: else:
projected = cv2.projectPoints(f.contour, f.rotation, f.translation, tracker.camera, tracker.dist_coeffs) projected = cv2.projectPoints(f.contour, f.rotation, f.translation, tracker.camera, tracker.dist_coeffs)
for [(x,y)] in projected[0]: for [(x,y)] in projected[0]:
@ -204,6 +210,14 @@ try:
packet.extend(bytearray(struct.pack("f", -z))) packet.extend(bytearray(struct.pack("f", -z)))
if not log is None: if not log is None:
log.write(f",{x},{-y},{-z}") log.write(f",{x},{-y},{-z}")
if f.current_features is None:
f.current_features = {}
for feature in features:
if not feature in f.current_features:
f.current_features[feature] = 0
packet.extend(bytearray(struct.pack("f", f.current_features[feature])))
if not log is None:
log.write(f",{f.current_features[feature]}")
if not log is None: if not log is None:
log.write("\r\n") log.write("\r\n")
log.flush() log.flush()

82
remedian.py Normal file
View file

@ -0,0 +1,82 @@
"""
Watch over a stream of numbers, incrementally learning their median.
Implemented via nested lists. New numbers are added to `lst[i]` and
when it fills up, it posts its median to `lst[i+1]`. Wen `lst[i+1]`
fills up, it posts the medians of its medians to `lst[i+2]`. Etc.
When a remedian is queried for the current median, it returns the
median of the last list with any numbers.
This approach is quite space efficient . E.g. four nested lists,
each with 64 items, require memory for 4*64 items yet can hold the
median of the median of the median of the median of over 17 million
numbers.
Example usage:
z=remedian()
for i in range(1000):
z + i
if not i % 100:
print(i, z.median())
Based on [The Remedian:A Robust Averaging Method for Large Data
Sets](http://web.ipac.caltech.edu/staff/fmasci/home/astro_refs/Remedian.pdf).
by Peter J. Rousseeuw and Gilbert W. Bassett Jr. Journal of the
American Statistical Association March 1990, Vol. 85, No. 409,
Theory and Methods
The code [remedianeg.py](remedianeg.py) compares this rig to just
using Python's built-in sort then reporing the middle number.
Assuming lists of length 64 and use of pypy3:
- Remedian is getting nearly as fast (within 20%) as raw sort after 500 items;
- While at the same time, avoids having to store all the numbers in RAM;
- Further, remedian's computed median is within 1% (or less) of the medians found via Python's sort.
_____
## Programmer's Guide
"""
# If `ordered` is `False`, do not sort `lst`
def median(lst,ordered=False):
assert lst,"median needs a non-empty list"
n = len(lst)
p = q = n//2
if n < 3:
p,q = 0, n-1
else:
lst = lst if ordered else sorted(lst)
if not n % 2: # for even-length lists, use mean of mid 2 nums
q = p -1
return lst[p] if p==q else (lst[p]+lst[q])/2
class remedian:
# Initialization
def __init__(i,inits=[], k=64, # after some experimentation, 64 works ok
about = None):
i.all,i.k = [],k
i.more,i._median=None,None
[i + x for x in inits]
# When full, push the median of current values to next list, then reset.
def __add__(i,x):
i._median = None
i.all.append(x)
if len(i.all) == i.k:
i.more = i.more or remedian(k=i.k)
i.more + i._medianPrim(i.all)
i.all = [] # reset
# If there is a next list, ask its median. Else, work it out locally.
def median(i):
return i.more.median() if i.more else i._medianPrim(i.all)
# Only recompute median if we do not know it already.
def _medianPrim(i,all):
if i._median == None:
i._median = median(all,ordered=False)
return i._median

View file

@ -8,6 +8,7 @@ import queue
import threading import threading
import copy import copy
from retinaface import RetinaFaceDetector from retinaface import RetinaFaceDetector
from remedian import remedian
def resolve(name): def resolve(name):
f = os.path.join(os.path.dirname(__file__), name) f = os.path.join(os.path.dirname(__file__), name)
@ -26,13 +27,13 @@ def clamp_to_im(pt, w, h):
y = h-1 y = h-1
return (int(x), int(y+1)) return (int(x), int(y+1))
def rotate(origin, point, angle): def rotate(origin, point, a):
angle = -angle a = -a
ox, oy = origin ox, oy = origin
px, py = point px, py = point
qx = ox + math.cos(angle) * (px - ox) - math.sin(angle) * (py - oy) qx = ox + math.cos(a) * (px - ox) - math.sin(a) * (py - oy)
qy = oy + math.sin(angle) * (px - ox) + math.cos(angle) * (py - oy) qy = oy + math.sin(a) * (px - ox) + math.cos(a) * (py - oy)
return qx, qy return qx, qy
def angle(p1, p2): def angle(p1, p2):
@ -45,10 +46,10 @@ def compensate(p1, p2):
a = angle(p1, p2) a = angle(p1, p2)
return rotate(p1, p2, a), a return rotate(p1, p2, a), a
def rotate_image(image, angle, center): def rotate_image(image, a, center):
(h, w) = image.shape[:2] (h, w) = image.shape[:2]
angle = np.rad2deg(angle) a = np.rad2deg(a)
M = cv2.getRotationMatrix2D((center[0], center[1]), angle, 1.0) M = cv2.getRotationMatrix2D((center[0], center[1]), a, 1.0)
rotated = cv2.warpAffine(image, M, (w, h)) rotated = cv2.warpAffine(image, M, (w, h))
return rotated return rotated
@ -115,6 +116,149 @@ def worker_thread(session, input, crop_info, queue, input_name, idx, tracker):
conf, lms = tracker.landmarks(output[0], crop_info) conf, lms = tracker.landmarks(output[0], crop_info)
queue.put((session, conf, lms, crop_info, idx)) queue.put((session, conf, lms, crop_info, idx))
class Feature():
def __init__(self, threshold=0.15, alpha=0.2, hard_factor=0.15, decay=0.001):
self.median = remedian()
self.min = None
self.max = None
self.hard_min = None
self.hard_max = None
self.threshold = threshold
self.alpha = alpha
self.hard_factor = hard_factor
self.decay = decay
self.last = 0
def update(self, x):
new = self.update_state(x)
filtered = self.last * self.alpha + new * (1 - self.alpha)
self.last = filtered
return filtered
def update_state(self, x):
self.median + x
median = self.median.median()
if self.min is None:
if x < median and (median - x) / median > self.threshold:
self.min = x
self.hard_min = self.min + self.hard_factor * (median - self.min)
return -1
return 0
else:
if x < self.min:
self.min = x
self.hard_min = self.min + self.hard_factor * (median - self.min)
return -1
if self.max is None:
if x > median and (x - median) / median > self.threshold:
self.max = x
self.hard_max = self.max - self.hard_factor * (self.max - median)
return 1
return 0
else:
if x > self.max:
self.max = x
self.hard_max = self.max - self.hard_factor * (self.max - median)
return 1
if self.min < self.hard_min:
self.min = self.hard_min * self.decay + self.min * (1 - self.decay)
if self.max > self.hard_max:
self.max = self.hard_max * self.decay + self.max * (1 - self.decay)
if x < median:
return - (1 - (x - self.min) / (median - self.min))
elif x > median:
return (x - median) / (self.max - median)
return 0
class FeatureExtractor():
def __init__(self):
self.eye_l = Feature()
self.eye_r = Feature()
self.eyebrow_updown_l = Feature()
self.eyebrow_updown_r = Feature()
self.eyebrow_quirk_l = Feature(threshold=0.05)
self.eyebrow_quirk_r = Feature(threshold=0.05)
self.eyebrow_steepness_l = Feature(threshold=0.05)
self.eyebrow_steepness_r = Feature(threshold=0.05)
self.mouth_corner_updown_l = Feature()
self.mouth_corner_updown_r = Feature()
self.mouth_corner_inout_l = Feature(threshold=0.02)
self.mouth_corner_inout_r = Feature(threshold=0.02)
self.mouth_open = Feature()
self.mouth_wide = Feature(threshold=0.02)
def align_points(self, a, b, pts):
a = tuple(a)
b = tuple(b)
alpha = angle(a, b)
alpha = np.rad2deg(alpha)
if alpha >= 90:
alpha = - (alpha - 180)
if alpha <= -90:
alpha = - (alpha + 180)
alpha = np.deg2rad(alpha)
aligned_pts = []
for pt in pts:
aligned_pts.append(np.array(rotate(a, pt, alpha)))
return alpha, np.array(aligned_pts)
def update(self, pts):
features = {}
norm_distance_x = np.mean([pts[0, 0] - pts[16, 0], pts[1, 0] - pts[15, 0]])
norm_distance_y = np.mean([pts[27, 1] - pts[28, 1], pts[28, 1] - pts[29, 1], pts[29, 1] - pts[30, 1]])
a1, f_pts = self.align_points(pts[42], pts[45], pts[[43, 44, 47, 46]])
f = np.clip((np.mean([f_pts[0,1], f_pts[1,1]]) - np.mean([f_pts[2,1], f_pts[3,1]])) / norm_distance_y, 0, None)
features["eye_l"] = self.eye_l.update(f)
a2, f_pts = self.align_points(pts[36], pts[39], pts[[37, 38, 41, 40]])
f = np.clip((np.mean([f_pts[0,1], f_pts[1,1]]) - np.mean([f_pts[2,1], f_pts[3,1]])) / norm_distance_y, 0, None)
features["eye_r"] = self.eye_r.update(f)
a3, _ = self.align_points(pts[0], pts[16], [])
a4, _ = self.align_points(pts[31], pts[35], [])
norm_angle = np.mean(list(map(np.rad2deg, [a1, a2, a3, a4])))
a, f_pts = self.align_points(pts[22], pts[26], pts[[22, 23, 24, 25, 26]])
features["eyebrow_steepness_l"] = self.eyebrow_steepness_l.update(-np.rad2deg(a) - norm_angle)
f = (np.mean([pts[22, 1], pts[26, 1]]) - pts[27, 1]) / norm_distance_y
features["eyebrow_updown_l"] = self.eyebrow_updown_l.update(f)
f = np.max(np.abs(np.array(f_pts[1:4]) - f_pts[0, 1])) / norm_distance_y
features["eyebrow_quirk_l"] = self.eyebrow_quirk_l.update(f)
a, f_pts = self.align_points(pts[17], pts[21], pts[[17, 18, 19, 20, 21]])
features["eyebrow_steepness_r"] = self.eyebrow_steepness_r.update(np.rad2deg(a) - norm_angle)
f = (np.mean([pts[17, 1], pts[21, 1]]) - pts[27, 1]) / norm_distance_y
features["eyebrow_updown_r"] = self.eyebrow_updown_r.update(f)
f = np.max(np.abs(np.array(f_pts[1:4]) - f_pts[0, 1])) / norm_distance_y
features["eyebrow_quirk_r"] = self.eyebrow_quirk_r.update(f)
upper_mouth_line = np.mean([pts[49, 1], pts[50, 1], pts[51, 1]])
center_line = np.mean([pts[50, 0], pts[60, 0], pts[27, 0], pts[30, 0], pts[64, 0], pts[55, 0]])
f = (upper_mouth_line - pts[62, 1]) / norm_distance_y
features["mouth_corner_updown_l"] = self.mouth_corner_updown_l.update(f)
f = abs(center_line - pts[62, 0]) / norm_distance_x
features["mouth_corner_inout_l"] = self.mouth_corner_inout_l.update(f)
f = (upper_mouth_line - pts[58, 1]) / norm_distance_y
features["mouth_corner_updown_r"] = self.mouth_corner_updown_r.update(f)
f = abs(center_line - pts[58, 0]) / norm_distance_x
features["mouth_corner_inout_r"] = self.mouth_corner_inout_r.update(f)
f = (np.mean([pts[59, 1], pts[60, 1], pts[61, 1]]) - np.mean([pts[65, 1], pts[64, 1], pts[63, 1]])) / norm_distance_y
features["mouth_open"] = self.mouth_open.update(f)
f = abs(pts[58, 0] - pts[62, 0]) / norm_distance_x
features["mouth_wide"] = self.mouth_wide.update(f)
return features
class FaceInfo(): class FaceInfo():
def __init__(self, id, tracker): def __init__(self, id, tracker):
self.id = id self.id = id
@ -128,7 +272,7 @@ class FaceInfo():
self.limit_3d_adjustment = True self.limit_3d_adjustment = True
self.update_count_delta = 75. self.update_count_delta = 75.
self.update_count_max = 2500. self.update_count_max = 7500.
def reset(self): def reset(self):
self.alive = False self.alive = False
@ -144,6 +288,9 @@ class FaceInfo():
self.pts_3d = None self.pts_3d = None
self.eye_blink = None self.eye_blink = None
self.bbox = None self.bbox = None
self.pnp_error = 0
self.features = FeatureExtractor()
self.current_features = {}
self.contour = np.zeros((21,3)) self.contour = np.zeros((21,3))
self.update_counts = np.zeros((66,2)) self.update_counts = np.zeros((66,2))
self.update_contour() self.update_contour()
@ -161,16 +308,22 @@ class FaceInfo():
self.contour = np.array(self.face_3d[self.contour_pts]) self.contour = np.array(self.face_3d[self.contour_pts])
def adjust_3d(self): def adjust_3d(self):
if self.conf < 0.4: if self.conf < 0.4 or self.pnp_error > 300:
return return
self.current_features = self.features.update(self.pts_3d[:, 0:2])
self.eye_blink = []
self.eye_blink.append(1 - min(max(0, -self.current_features["eye_r"]), 1))
self.eye_blink.append(1 - min(max(0, -self.current_features["eye_l"]), 1))
max_runs = 1 max_runs = 1
eligible = np.arange(0, 66) eligible = np.arange(0, 66)
changed_any = False changed_any = False
update_type = -1
d_o = np.ones((66,)) d_o = np.ones((66,))
d_c = np.ones((66,)) d_c = np.ones((66,))
for runs in range(max_runs): for runs in range(max_runs):
r = 1.0 + np.random.random_sample((66,3)) * 0.02 - 0.01 r = 1.0 + np.random.random_sample((66,3)) * 0.02 - 0.01
update_type = -1
if self.euler[0] > -165 and self.euler[0] < 145: if self.euler[0] > -165 and self.euler[0] < 145:
continue continue
elif self.euler[1] > -10 and self.euler[1] < 20: elif self.euler[1] > -10 and self.euler[1] < 20:
@ -226,7 +379,10 @@ class FaceInfo():
weights[:, :] = self.lms[0:66, 2:3] weights[:, :] = self.lms[0:66, 2:3]
weights[weights > 0.7] = 1.0 weights[weights > 0.7] = 1.0
weights = 1.0 - weights weights = 1.0 - weights
self.face_3d[0:66] = self.face_3d[0:66] * weights + updated[0:66] * (1. - weights) update_indices = np.arange(0, 66)
if self.limit_3d_adjustment:
update_indices = np.nonzero(self.update_counts[:, update_type] <= self.update_count_max)[0]
self.face_3d[update_indices] = self.face_3d[update_indices] * weights[update_indices] + updated[update_indices] * (1. - weights[update_indices])
self.update_contour() self.update_contour()
class Tracker(): class Tracker():
@ -391,6 +547,7 @@ class Tracker():
self.no_gaze = no_gaze self.no_gaze = no_gaze
self.debug_gaze = False self.debug_gaze = False
self.face_info = [FaceInfo(id, self) for id in range(max_faces)] self.face_info = [FaceInfo(id, self) for id in range(max_faces)]
self.fail_count = 0
def landmarks(self, tensor, crop_info): def landmarks(self, tensor, crop_info):
crop_x1, crop_y1, scale_x, scale_y = crop_info crop_x1, crop_y1, scale_x, scale_y = crop_info
@ -484,12 +641,19 @@ class Tracker():
pts_3d[i,:] = pt_3d[:] pts_3d[i,:] = pt_3d[:]
pnp_error = np.sqrt(pnp_error / (2.0 * image_pts.shape[0])) pnp_error = np.sqrt(pnp_error / (2.0 * image_pts.shape[0]))
if pnp_error > 50: if pnp_error > 300:
# Something went wrong with adjusting the 3D model self.fail_count += 1
face_info.face_3d = self.face_3d if self.fail_count > 5:
face_info.rotation = np.array([0.0, 0.0, 0.0], np.float32) # Something went wrong with adjusting the 3D model
face_info.translation = np.array([0.0, 0.0, 0.0], np.float32) if not self.silent:
face_info.update_counts = np.zeros((66,2)) print(f"Detected anomaly when 3D fitting face {face_info.id}. Resetting.")
face_info.face_3d = copy.copy(self.face_3d)
face_info.rotation = None
face_info.translation = np.array([0.0, 0.0, 0.0], np.float32)
face_info.update_counts = np.zeros((66,2))
face_info.update_contour()
else:
self.fail_count = 0
euler = cv2.RQDecomp3x3(rmat)[0] euler = cv2.RQDecomp3x3(rmat)[0]
return True, matrix_to_quaternion(rmat), euler, pnp_error, pts_3d, lms return True, matrix_to_quaternion(rmat), euler, pnp_error, pts_3d, lms
@ -516,20 +680,20 @@ class Tracker():
((cx1, cy1), (cx2, cy2)) = corners ((cx1, cy1), (cx2, cy2)) = corners
c1 = np.array([cx1, cy1]) c1 = np.array([cx1, cy1])
c2 = np.array([cx2, cy2]) c2 = np.array([cx2, cy2])
c2, angle = compensate(c1, c2) c2, a = compensate(c1, c2)
center = (c1 + c2) / 2.0 center = (c1 + c2) / 2.0
radius = np.linalg.norm(c1 - c2) / 2.0 radius = np.linalg.norm(c1 - c2) / 2.0
radius = np.array([radius * 1.4, radius * 1.2]) radius = np.array([radius * 1.4, radius * 1.2])
upper_left = clamp_to_im(center - radius, w, h) upper_left = clamp_to_im(center - radius, w, h)
lower_right = clamp_to_im(center + radius, w, h) lower_right = clamp_to_im(center + radius, w, h)
return upper_left, lower_right, center, radius, c1, angle return upper_left, lower_right, center, radius, c1, a
def prepare_eye(self, frame, full_frame, lms, flip): def prepare_eye(self, frame, full_frame, lms, flip):
outer_pt = tuple(lms[0]) outer_pt = tuple(lms[0])
inner_pt = tuple(lms[1]) inner_pt = tuple(lms[1])
h, w, _ = frame.shape h, w, _ = frame.shape
(x1, y1), (x2, y2), center, radius, reference, angle = self.corners_to_eye((outer_pt, inner_pt), w, h, flip) (x1, y1), (x2, y2), center, radius, reference, a = self.corners_to_eye((outer_pt, inner_pt), w, h, flip)
im = rotate_image(frame[:, :, ::], angle, reference) im = rotate_image(frame[:, :, ::], a, reference)
im = im[int(y1):int(y2), int(x1):int(x2),:] im = im[int(y1):int(y2), int(x1):int(x2),:]
if np.prod(im.shape) < 1: if np.prod(im.shape) < 1:
return None, None, None, None, None, None return None, None, None, None, None, None
@ -546,7 +710,7 @@ class Tracker():
im = im.astype(np.float32)[:,:,::-1] / self.std - self.mean im = im.astype(np.float32)[:,:,::-1] / self.std - self.mean
im = np.expand_dims(im, 0) im = np.expand_dims(im, 0)
im = np.transpose(im, (0,3,2,1)) im = np.transpose(im, (0,3,2,1))
return im, x1, y1, scale, reference, angle return im, x1, y1, scale, reference, a
def extract_face(self, frame, lms): def extract_face(self, frame, lms):
lms = np.array(lms)[:,0:2][:,::-1] lms = np.array(lms)[:,0:2][:,::-1]
@ -572,10 +736,10 @@ class Tracker():
e_y = [0,0] e_y = [0,0]
scale = [0,0] scale = [0,0]
reference = [None, None] reference = [None, None]
angle = [0, 0] angles = [0, 0]
face_frame, lms, offset = self.extract_face(frame, lms) face_frame, lms, offset = self.extract_face(frame, lms)
(right_eye, e_x[0], e_y[0], scale[0], reference[0], angle[0]) = self.prepare_eye(face_frame, frame, np.array([lms[36,0:2], lms[39,0:2]]), False) (right_eye, e_x[0], e_y[0], scale[0], reference[0], angles[0]) = self.prepare_eye(face_frame, frame, np.array([lms[36,0:2], lms[39,0:2]]), False)
(left_eye, e_x[1], e_y[1], scale[1], reference[1], angle[1]) = self.prepare_eye(face_frame, frame, np.array([lms[42,0:2], lms[45,0:2]]), True) (left_eye, e_x[1], e_y[1], scale[1], reference[1], angles[1]) = self.prepare_eye(face_frame, frame, np.array([lms[42,0:2], lms[45,0:2]]), True)
if right_eye is None or left_eye is None: if right_eye is None or left_eye is None:
return [(1.0, 0.0, 0.0, 0.0), (1.0, 0.0, 0.0, 0.0)] return [(1.0, 0.0, 0.0, 0.0), (1.0, 0.0, 0.0, 0.0)]
both_eyes = np.concatenate((right_eye, left_eye)) both_eyes = np.concatenate((right_eye, left_eye))
@ -617,7 +781,7 @@ class Tracker():
else: else:
eye_x = e_x[i] + scale[i][0] * (32. - eye_x) eye_x = e_x[i] + scale[i][0] * (32. - eye_x)
eye_y = e_y[i] + scale[i][1] * eye_y eye_y = e_y[i] + scale[i][1] * eye_y
eye_x, eye_y = rotate(reference[i], (eye_x, eye_y), -angle[i]) eye_x, eye_y = rotate(reference[i], (eye_x, eye_y), -angles[i])
eye_x = eye_x + offset[0] eye_x = eye_x + offset[0]
eye_y = eye_y + offset[1] eye_y = eye_y + offset[1]
@ -794,7 +958,7 @@ class Tracker():
bbox = (min_y, min_x, max_y - min_y, max_x - min_x) bbox = (min_y, min_x, max_y - min_y, max_x - min_x)
detected.append(bbox) detected.append(bbox)
face_info.bbox = bbox face_info.bbox = bbox
face_info.eye_blink = (eye_state[0][0], eye_state[1][0]) #face_info.eye_blink = (eye_state[0][0], eye_state[1][0])
results.append(face_info) results.append(face_info)
if len(detected) > 0: if len(detected) > 0: