Gradient calculation woes

Hi all, ultimately I’m trying to simulate Chladni patterns (made by vibrating a metal plate with particles on it). I’ve come across this equation that gives Chladni patterns (a, b, m, n are parameters).

double ofApp::f(double a, double b, double m, double n, double x, double y){
	
	double zpos = a*glm::sin(3.14*x*n)*glm::sin(3.14*y*m) + b*glm::sin(3.14*x*m)*glm::sin(3.14*y*n);
	return zpos;
	
}

(small aside, things break when I try to use the PI constant defined in ofMathConstants.h. It also breaks if I manually include more digits like 3.14159).

I’m trying to take the gradient of this function, and use it to push particles around (into the nodes of the function f(a, b, m, n, x, y) ideally). I do this as below:

ofVec2f ofApp::gradCalc(double a, double b, double m, double n, double x, double y){
	
	ofVec2f v;
	double eps = 1e-10;
	double dx = eps;
	double dfdx = (f(a, b, m, n, x+dx, y) - f(a, b, m, n, x, y)) / (dx);
	double dy = eps;
	double dfdy = (f(a, b, m, n, x, y+dy) - f(a, b, m, n, x, y)) / (dy);
	
	v.set(dfdx, dfdy);
	return v;

}

I can see that this is working… almost. I can draw the gradient field and it looks how I’d expect. However when I try to use this field to move particles, I run into trouble. To debug I focused on just one particle, and I set it up so that I can click on the screen and it’ll show what the gradient is at that point. Sometimes this will look correct and align with the field, but other times it looks very wrong. The magnitude always looks correct though.

Here it is looking correct:

But at a position very close by I get unexpected behavior:

I was thinking it had to do with numerical precision so I went through and made all floats doubles. I’m not very versed in C++ so I feel like I’m overlooking something silly!

Thanks for reading :slight_smile:

Here’s the full myApp.cpp and header:

#include "ofApp.h"

vector<ofVec3f> pos;
vector<ofVec3f> origPos;
ofEasyCam cam;
ofNode node;
double a = 5;
double b = 5;
double m = 5;
double n = 5;
ofVec3f mousePos (0,0,0);

bool framestep = false;

//--------------------------------------------------------------
void ofApp::setup(){
	
	node.setPosition(ofGetWidth()/2, ofGetHeight()/2, 0);
	cam.setTarget(node);
	
	ofSetFrameRate( 60 );
		
	for(double i = 0; i < 50; i++){
		for(double j = 0; j < 50; j++){
			ofVec3f v;
			double zpos = f(a, b, m, n, i*20, j*20);
			//cout << zpos << endl;
			v.set(i*20, j*20, zpos);
		
			pos.push_back(v);
			
		}
	}
	
	origPos.assign(pos.begin(), pos.end());
	
}

//--------------------------------------------------------------
void ofApp::update(){
	
	if (framestep == true){
		
		cout << mousePos << endl;
		cout << gradCalc(a, b, m, n, mousePos.x, mousePos.y) << endl;
		
		framestep = false;
		
	}
	
	 
}

//--------------------------------------------------------------
void ofApp::draw(){
	
	//try to draw the gradient vectors to see if that calculation's working
	
	//cam.begin();
	
        //this loop draws the gradient field on the 50x50 grid
	for(int i = 0; i < pos.size(); i++){
		ofSetColor(0, 0, 0);
		ofDrawArrow(origPos[i], origPos[i] + gradCalc(a, b, m, n, origPos[i].x, origPos[i].y), 5);
	}
	
        //draw the gradient at the point of mouse click
	ofSetColor(255, 255, 255);
	ofDrawCircle(mousePos.x, mousePos.y, 5);
	ofDrawArrow(mousePos, mousePos + gradCalc(a, b, m, n, mousePos.x, mousePos.y), 5);
	
	//cam.end();
	
	string fpsStr = "frame rate: "+ofToString(ofGetFrameRate(), 2);
	ofDrawBitmapString(fpsStr, 100,100);

}

double ofApp::f(double a, double b, double m, double n, double x, double y){
	
	double zpos = a*glm::sin(3.14*x*n)*glm::sin(3.14*y*m) + b*glm::sin(3.14*x*m)*glm::sin(3.14*y*n);
	return zpos;
	
}

ofVec2f ofApp::gradCalc(double a, double b, double m, double n, double x, double y){
	
	ofVec2f v;
	double eps = 1e-10;
	double dx = eps;
	double dfdx = (f(a, b, m, n, x+dx, y) - f(a, b, m, n, x, y)) / (dx);
	double dy = eps;
	double dfdy = (f(a, b, m, n, x, y+dy) - f(a, b, m, n, x, y)) / (dy);
	
	v.set(dfdx, dfdy);
	return v;

}

//--------------------------------------------------------------
void ofApp::keyPressed(int key){

}

//--------------------------------------------------------------
void ofApp::keyReleased(int key){

}

//--------------------------------------------------------------
void ofApp::mouseMoved(int x, int y ){

}

//--------------------------------------------------------------
void ofApp::mouseDragged(int x, int y, int button){

}

//--------------------------------------------------------------
void ofApp::mousePressed(int x, int y, int button){
	
}

//--------------------------------------------------------------
void ofApp::mouseReleased(int x, int y, int button){
	
	mousePos.set(ofGetMouseX(), ofGetMouseY());
	framestep = true;

}

//--------------------------------------------------------------
void ofApp::mouseEntered(int x, int y){

}

//--------------------------------------------------------------
void ofApp::mouseExited(int x, int y){

}

//--------------------------------------------------------------
void ofApp::windowResized(int w, int h){

}

//--------------------------------------------------------------
void ofApp::gotMessage(ofMessage msg){

}

//--------------------------------------------------------------
void ofApp::dragEvent(ofDragInfo dragInfo){ 

}

header:

#pragma once

#include "ofMain.h"

class ofApp : public ofBaseApp{

	public:
		void setup();
		void update();
		void draw();

		void keyPressed(int key);
		void keyReleased(int key);
		void mouseMoved(int x, int y );
		void mouseDragged(int x, int y, int button);
		void mousePressed(int x, int y, int button);
		void mouseReleased(int x, int y, int button);
		void mouseEntered(int x, int y);
		void mouseExited(int x, int y);
		void windowResized(int w, int h);
		void dragEvent(ofDragInfo dragInfo);
		void gotMessage(ofMessage msg);
	
		ofVec2f gradCalc(double a, double b, double m, double n, double x, double y);
		double f(double a, double b, double m, double n, double x, double y);
		
};

Try changing the last line of your function gradCalc to:

return v*0.1;

I think the issue is your eye is drawn to the arrow heads, but if you follow the base of the arrow which is actually the location being calculated it’s quite a different spot.

So in actuality you are right on a seam between the two tiles, you just can’t tell because the arrows are so long. ( I might be wrong here ).

I’ve once had an issue multiplying recursively small numbers to make an “easing” function, and when the number was near zero I had some unexpected behavior.

I’ve solving this substituting the old recursive division for ofLerp, to interpolate old and new values, something like this:

valEasy = ofLerp(valEasy, val, easing);

PS: Not sure if it is the same issue you are experiencing

hmm I thought this was right at first, but I think there’s something else going on. In this clip I’m going straight up hill, but the direction still flip-flops

I think it is related to small numbers or precision or something, but honestly idk. I’ve shown the limit of my numerics knowledge already haha. lerping or truncating something may help

not a “solution” per se but some ideas:

interesting that 3.141592 breaks things, this is the reddest flag. I think your 3.14 is not really in a degree context and behaves more like an arbitrary multiplier which happens to look nice at 3.14, but it is very sensitive to little movements – perhaps and what we are seeing is actually a form of aliasing of a different frequency and the scale is super off.

note that ofVec3f is glm::vec3 and glm::vec* are floats and operations such as glm::sin are returning floats. all this to say using double will not affect the behaviour. for sure all your initial positions (both mouseclicks and as generated from the code) are integer values.

also you can retrieve your platforms actual float epsilon for a type like this: std::numeric_limits<float>::epsilon();

I would recommend making a GUI with your params and playing with them to understand better the relationship of units with the visual result.

1 Like

you’re right on, the natural scale for the function as I’ve written above is like [-0.5, 0.5], not 0 to 1000… the aliasing really threw me for a loop, and I’m still really surprised my original grid spacing aligned so neatly! thanks for your comments. this is what I get for not plotting it first