Here is my second take using your approach rather than the one i described earlier.
Before going into the process of randomizing the shape, I started to clean a few things:

your angle step (0.0001) is really small so first it takes way more time to produce the shape and second it can actually produce some artifacts since you are asking processing to draw lines between points that are on the same pixel. I lowered it to 0.01 and it should be more than enough.

I changed the way you calculate the (x, y) coordinate of your shapes because it is quite confusing. It’s almost polar coordinates but it’s not. And because of that the meaning of the angle is not intuitive (I lost quite some times figuring this out actually).
The way I did it was to first log all the (x, y) coordinates of one of the 4 shapes. Based on those values, i calculated the angle and the radius of the point in polar coordinates. The next part was plotting the (theta, r) coordinates and find a function that nicely approximate the shape. I did it in Excel and it looks like this for one of the shape:
The blue is the original and the orange the fitting function.
In the code, this is the purpose of the functions f
and g
With this fresh start, the next thing to do was to produce a smooth noise with a period of 2 PI since the idea is to slightly change the r value for each angles. This is what the PRNG
class is doing.
Finally, I don’t want to change the r value on the top part and the bottom part but only on the side. The way around this is to create a masking function that will be 0 for the angle corresponding to the bottom and top parts and 1 for the angle corresponding to the sides. The is the purpose of the scaleFactor
function. The function itself look like this:
Here’s the result:
Code
PRNG prng = new PRNG(6, TWO_PI);
void setup() {
size(800, 700, P2D);
background(20);
for (int c = 0; c < width  100; c = c + 120) {
for (int d = 0; d < 300; d = d + 290) {
vagina(c + 100, d + 200);
}
}
}
float scaleFactor(float x, float a, float b) {
x = ((x + HALF_PI) % PI)  HALF_PI;
float c = 0;
return 1 / (1 + pow(abs((x c) / a), 2 * b));
}
float f(float x, float amp, float power, float offset) {
return amp * (2  pow(abs(cos(x)), power)) + offset;
}
float g(float x, float scale, float amp, float xOffset, float yOffset) {
x = abs(((x + HALF_PI) % PI)  HALF_PI);
return (amp / (log10(scale * x +xOffset))) + yOffset;
}
float log10 (float x) {
return (log(x) / log(10));
}
void vagina(float c, float d) {
stroke(230);
strokeWeight(1);
noFill();
prng.randomSeed(millis());
beginShape();
for (float a = 0; a < TWO_PI; a = a + 0.01) {
float r = f(a, 50, 0.4, 0);
r += (prng.smoothNoise(a)  0.5) * 8 * scaleFactor(a, HALF_PI / 2.0, 6);
float x = c + r * cos(a);
float y = d + r * sin(a);
vertex(x, y);
}
endShape();
prng.randomSeed(millis());
beginShape();
for (float a = 0; a < TWO_PI; a = a + 0.01) {
float r = f(a, 100, 0.28, 55);
r += (prng.smoothNoise(a)  0.5) * 10 * scaleFactor(a, HALF_PI / 2.0, 6);
float x = c + r * cos(a);
float y = d + r * sin(a);
vertex(x, y);
}
endShape();
prng.randomSeed(millis());
beginShape();
for (float a = 0; a < TWO_PI; a = a + 0.01) {
float r = f(a, 120, 0.25, 100);
r += (prng.smoothNoise(a)  0.5) * 4 * scaleFactor(a, HALF_PI / 2.0, 6);
float x = c + r * cos(a);
float y = d + r * sin(a);
vertex(x, y);
}
endShape();
prng.randomSeed(millis());
beginShape();
for (float a = 0; a < TWO_PI; a = a + 0.01) {
float r = g(a, 2.48, 6, 5, 3.5);
r += (prng.smoothNoise(a)  0.5) * 2 * scaleFactor(a, HALF_PI / 2.0, 6);
float x = c + r * cos(a);
float y = d + r * sin(a);
vertex(x, y);
}
endShape();
}
class PRNG {
private long seed;
private long a;
private long c;
private long m32;
private float[] samplePts;
private int detail;
private float period;
PRNG(int detail, float period) {
this(System.currentTimeMillis(), detail, period);
}
PRNG(long seed, int detail, float period) {
this.a = 1664525;
this.c = 1013904223;
this.m32 = 0xFFFFFFFFL;
this.detail = detail;
this.period = period;
this.randomSeed(seed);
}
void randomSeed(long newSeed) {
this.seed = newSeed%Integer.MAX_VALUE;
initSmoothNoise();
}
long nextLong() {
this.seed = this.seed * a + c & m32;
return this.seed;
}
int nextInt() {
return (int)(this.nextLong()%Integer.MAX_VALUE);
}
float random() {
return random(0, 1);
}
float random(float max) {
return random(0, max);
}
float random(float min, float max) {
return map(this.nextInt(), 0, Integer.MAX_VALUE, min, max);
}
private void initSmoothNoise() {
samplePts = new float[detail + 1];
for (int i = 0; i < samplePts.length; i++) {
samplePts[i] = random();
}
}
private float lerpValues(float a, float b, float t) {
return a + fade(t) * (b  a);
}
private float fade(float t) {
return t * t * t * ( t * ( t * 6  15 ) + 10 );
}
void changePeriod(float p) {
period = p;
initSmoothNoise();
}
void changeDetail(int d) {
detail = d;
initSmoothNoise();
}
float smoothNoise(float x) {
float i = map((x % period), 0, period, 0, detail + 1);
int prev = floor(i);
int next = (int)((i + 1) % (detail + 1));
float t = i  prev;
return lerpValues(samplePts[prev], samplePts[next], t);
}
}