处理方程优化反应扩散java

Question

我正在用java（和Processing）处理Gray Scott patterns 我想优化它。 我没有任何GPU加速技能。因此，我尝试优化方程：

U[i][j] += ( U*lapU - uvv + F*(1 - u) ) * 1.38 ;
V[i][j] += ( V*lapV + uvv - (K+F)*v   ) * 0.63 ;

但也许还有其他的方法，这里是函数

如果你有任何提高速度的想法

float uvv, u, v, diffU, diffV, F, K, lapU, lapV;
float[] MINI = { 0.00, 0.01, 0.03, 0.005 }; 
float[] MAXI = { 0.15, 0.08, 0.11, 0.05 }; 

PImage reactionDiffusion (PImage img) {

  int W = img.width, H = img.height;  
  int[][] offsetW = new int[W][2], offsetH = new int[H][2];
  float[][]  U = new float[W][H],  V = new float[W][H];
  float time = millis();

  //  INITIALISATION
  for (int i = 0; i < W; ++i) {
    for (int j = 0; j < H; ++j) {
      U[i][j] = 0.15 ;
      V[i][j] = 0.7  ;
    }
  }  

  // init parameters
  float[][][] fkuv = new float[W][H][4];  

  for (int i = 0; i<W; ++i){
    for (int j = 0; j<H; ++j){
        fkuv[i][j][0] = map( i, 0, H, MINI[0], MAXI[0] );
        fkuv[i][j][1] = map( j, 0, W, MAXI[1], MINI[0] );  
        fkuv[i][j][2] = map(150,0,200,0,MAXI[2]);
        fkuv[i][j][3] = map(50,0,200,0,MAXI[3]);
    }
  }
  //Set up offsets
  for (int i=0; i < W; ++i) { offsetW[i][0] = i-1; offsetW[i][1] = i+1; }
  for (int i=0; i < H; ++i) { offsetH[i][0] = i-1; offsetH[i][1] = i+1; }
  offsetW[0][0] = 0; offsetW[W-1][1] = W-1;
  offsetH[0][0] = 0; offsetH[H-1][1] = H-1;

  for (int n = 0; n< 2000 ; ++n){ // itérations
    for (int i = 0; i < W; ++i) {
      for (int j = 0; j < H; ++j) {

        F = fkuv[i][j][0] ;
        K = fkuv[i][j][1] ;

        u = U[i][j];  
        v = V[i][j]; 

        uvv = u*v*v;
        lapU = U[offsetW[i][0]][j] +U[offsetW[i][1]][j] +U[i][offsetH[j][0]] +U[i][offsetH[j][1]] -4*u;
        lapV = V[offsetW[i][0]][j] +V[offsetW[i][1]][j] +V[i][offsetH[j][0]] +V[i][offsetH[j][1]] -4*v;

        U[i][j] += ( U*lapU - uvv + F*(1 - u) ) * 1.38 ;
        V[i][j] += ( V*lapV + uvv - (K+F)*v   ) * 0.63 ;
      } 
    }
  }
  img.loadPixels();
    int pShift;
    for (int i = 0; i < W; i++) {
      for (int j = 0; j < H; j++) {
        pShift = int( U[i][j]*255 ) ;
        img.pixels[j*W+i] = 0xff000000 | (pShift << 16) | (pShift << 8) | pShift  ;
      }
    }
  img.updatePixels();

  println("time : "+ ( millis()-time ) /1000) ; 

  return img;
}

谢谢

Answer 1

您可以在this answer中找到一些优化技巧

在GPU加速方面，你不需要学习所有东西。您应该能够理解顶点着色器和片段着色器之间的区别以及如何使用统一变量。看看这个awesome PShader tutorial

处理中也有片段PShader示例。 一旦掌握了窍门，就可以将反应扩散从shadertoy移植到处理中（只需非常小的调整）

下面是一个非常粗略的移植this shader示例

//ReactionDiffusion
// https://www.shadertoy.com/view/MdKXzw

PShader bufA,shader;

void setup(){
  size(640,480,P2D);
  noStroke();

  //uniform sampler2D iChannel0;
  //uniform vec2 resolution;
  //uniform vec2 mouse;
  //uniform float time;
  bufA = loadShader("BufA.frag");
  bufA.set("resolution",(float)width,(float)height);
  //uniform sampler2D iChannel0;
  //uniform vec3 iChannelResolution[4];
  //uniform vec2 resolution;
  //uniform float time;
  shader = loadShader("shader.frag");
  shader.set("resolution",(float)width,(float)height);
  shader.set("iChannelResolution[0]",(float)width,(float)height,72);

}
void draw(){

  bufA.set("time",frameCount * .1);
  bufA.set("mouse",(float)mouseX,(float)mouseY,mousePressed ? 1 : 0,mousePressed && mouseButton == RIGHT ? 1 : 0);

  shader(bufA);
  background(0);
  rect(0,0,width,height);
  bufA.set("iChannel0",get());
  //2nd pass
  shader.set("iChannel0",get());
  shader.set("time",frameCount * .1);
  shader(shader);
  rect(0,0,width,height);
}

BufA。frag

uniform sampler2D iChannel0;
uniform vec2 resolution;
uniform vec4 mouse;
uniform float time;

vec2 vv(vec2 uv, vec2 duv, vec2 texel) { 
    return texture2D(iChannel0, uv + (duv * texel)).xy; 
}

vec2 laplacian(vec2 ab, vec2 uv, vec2 texel) {
    return 
        ab * -1. +
        vv(uv, vec2(0., -1.), texel) * .2 +
        vv(uv, vec2(0., 1.), texel) * .2 +
        vv(uv, vec2(1., 0.), texel) * .2 +
        vv(uv, vec2(-1., 0.), texel) * .2 +
        vv(uv, vec2(-1., -1.), texel) * .05 +
        vv(uv, vec2(-1., 1.), texel) * .05 +
        vv(uv, vec2(1., -1.), texel) * .05 +
        vv(uv, vec2(1., 1.), texel) * .05;
}

void main()
{

    vec2 uv = gl_FragCoord.xy / resolution.xy;
    vec2 texel = 1. / resolution.xy;
    vec4 m = mouse / resolution.xyxy;

    vec2 c = vec2(.5+cos(time)*0.3,.5+sin(time)*0.2);
    float lc = length(uv-c);

    if (m.z > 0. && length(m.xy-uv) < .02) {
        gl_FragColor = vec4(0., 1., 0., 1.);
    } else if (lc < .02-time*.003) {
        gl_FragColor = vec4(0., 1., 0., 1.);
    } else {

        vec2 ab = vv(uv, vec2(0., 0.), texel);
        vec2 lp = laplacian(ab, uv, texel);
        float a = ab.x;
        float b = ab.y;
        float la = lp.x;
        float lb = lp.y;
        float f = .028+(uv.y+c.y)*.0035;
        float k = .055+(uv.x+c.x)*.003;
        float dx = mod((time+uv.x+uv.y)*.5, .8);
        float da = .4+dx;
        float db = .2+dx*.5;


        float na = a + (da * la - a*b*b + f*(1.-a));
        float nb = b + (db * lb + a*b*b - (k+f)*b);

        gl_FragColor = vec4(na, nb, 0., 1.);
    }

}

着色器。frag

uniform sampler2D iChannel0;
uniform vec3 iChannelResolution[4];
uniform vec2 resolution;
uniform float time;

struct C_Sample
{
    vec3 vAlbedo;
    vec3 vNormal;
};

C_Sample SampleMaterial(const in vec2 vUV, sampler2D sampler,  const in vec2 vTextureSize, const in float fNormalScale)
{
    C_Sample result;

    vec2 vInvTextureSize = vec2(1.0) / vTextureSize;

    vec3 cSampleNegXNegY = clamp(texture2D(sampler, vUV + (vec2(-1.0, -1.0)) * vInvTextureSize.xy).ggb*2., 0., 1.);
    vec3 cSampleZerXNegY = clamp(texture2D(sampler, vUV + (vec2( 0.0, -1.0)) * vInvTextureSize.xy).ggb*2., 0., 1.);
    vec3 cSamplePosXNegY = clamp(texture2D(sampler, vUV + (vec2( 1.0, -1.0)) * vInvTextureSize.xy).ggb*2., 0., 1.);

    vec3 cSampleNegXZerY = clamp(texture2D(sampler, vUV + (vec2(-1.0, 0.0)) * vInvTextureSize.xy).ggb*2., 0., 1.);
    vec3 cSampleZerXZerY = clamp(texture2D(sampler, vUV + (vec2( 0.0, 0.0)) * vInvTextureSize.xy).ggb*2., 0., 1.);
    vec3 cSamplePosXZerY = clamp(texture2D(sampler, vUV + (vec2( 1.0, 0.0)) * vInvTextureSize.xy).ggb*2., 0., 1.);

    vec3 cSampleNegXPosY = clamp(texture2D(sampler, vUV + (vec2(-1.0,  1.0)) * vInvTextureSize.xy).ggb*2., 0., 1.);
    vec3 cSampleZerXPosY = clamp(texture2D(sampler, vUV + (vec2( 0.0,  1.0)) * vInvTextureSize.xy).ggb*2., 0., 1.);
    vec3 cSamplePosXPosY = clamp(texture2D(sampler, vUV + (vec2( 1.0,  1.0)) * vInvTextureSize.xy).ggb*2., 0., 1.);

    // convert to linear    
    vec3 cLSampleNegXNegY = cSampleNegXNegY * cSampleNegXNegY;
    vec3 cLSampleZerXNegY = cSampleZerXNegY * cSampleZerXNegY;
    vec3 cLSamplePosXNegY = cSamplePosXNegY * cSamplePosXNegY;

    vec3 cLSampleNegXZerY = cSampleNegXZerY * cSampleNegXZerY;
    vec3 cLSampleZerXZerY = cSampleZerXZerY * cSampleZerXZerY;
    vec3 cLSamplePosXZerY = cSamplePosXZerY * cSamplePosXZerY;

    vec3 cLSampleNegXPosY = cSampleNegXPosY * cSampleNegXPosY;
    vec3 cLSampleZerXPosY = cSampleZerXPosY * cSampleZerXPosY;
    vec3 cLSamplePosXPosY = cSamplePosXPosY * cSamplePosXPosY;

    result.vAlbedo = cLSampleNegXNegY;  

    vec3 vScale = vec3(0.3333);

    float fSampleNegXNegY = dot(cLSampleNegXNegY, vScale);
    float fSampleZerXNegY = dot(cLSampleZerXNegY, vScale);
    float fSamplePosXNegY = dot(cLSamplePosXNegY, vScale);

    float fSampleNegXZerY = dot(cLSampleNegXZerY, vScale);
    float fSampleZerXZerY = dot(cLSampleZerXZerY, vScale);
    float fSamplePosXZerY = dot(cLSamplePosXZerY, vScale);

    float fSampleNegXPosY = dot(cLSampleNegXPosY, vScale);
    float fSampleZerXPosY = dot(cLSampleZerXPosY, vScale);
    float fSamplePosXPosY = dot(cLSamplePosXPosY, vScale);


    // Sobel operator - http://en.wikipedia.org/wiki/Sobel_operator

    vec2 vEdge;
    vEdge.x = (fSampleNegXNegY - fSamplePosXNegY) * 0.25 
            + (fSampleNegXZerY - fSamplePosXZerY) * 0.5
            + (fSampleNegXPosY - fSamplePosXPosY) * 0.25;

    vEdge.y = (fSampleNegXNegY - fSampleNegXPosY) * 0.25 
            + (fSampleZerXNegY - fSampleZerXPosY) * 0.5
            + (fSamplePosXNegY - fSamplePosXPosY) * 0.25;

    result.vNormal = normalize(vec3(vEdge * fNormalScale, 1.0));    

    return result;
}

void main()
{   
    vec2 vUV = gl_FragCoord.xy / resolution.xy;

    C_Sample materialSample;

    float fNormalScale = 10.0;
    materialSample = SampleMaterial( vUV, iChannel0, iChannelResolution[0].xy, fNormalScale );

    float fLightHeight = 0.5;
    float fViewHeight = 2.0;

    vec3 vSurfacePos = vec3(vUV, 0.0);

    vec3 vViewPos = vec3(0.5, 0.5, fViewHeight);

    vec3 vLightPos = vec3( vec2(sin(time*.5),cos(time*.5)) * 0.25 + 0.5 , fLightHeight);

    vec3 vDirToView = normalize( vViewPos - vSurfacePos );
    vec3 vDirToLight = normalize( vLightPos - vSurfacePos );

    float fNDotL = clamp( dot(materialSample.vNormal, vDirToLight), 0.0, 1.0);
    float fDiffuse = fNDotL;

    vec3 vHalf = normalize( vDirToView + vDirToLight );
    float fNDotH = clamp( dot(materialSample.vNormal, vHalf), 0.0, 1.0);
    float fSpec = pow(fNDotH, 10.0) * fNDotL * 0.5;

    vec3 vResult = materialSample.vAlbedo * fDiffuse + fSpec;

    vResult = sqrt(vResult);

    gl_FragColor = vec4(vResult,1.0);
}

此外，您应该签出toxiclibs。它已经实现了与simutils包中的GrayScott类相同的文件

下面是一个例子：

/**
 * <p>GrayScottToneMap shows how to use the ColorGradient & ToneMap classes of the
 * colorutils package to create a tone map for rendering the results of
 * the Gray-Scott reaction-diffusion.</p>
 *
 * <p><strong>Usage:</strong><ul>
 * <li>click + drag mouse to draw dots used as simulation seed</li>
 * <li>press any key to reset</li>
 * </ul></p>
 */

/* 
 * Copyright (c) 2010 Karsten Schmidt
 * 
 * This demo & library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * 
 * http://creativecommons.org/licenses/LGPL/2.1/
 * 
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

import toxi.sim.grayscott.*;
import toxi.math.*;

import toxi.color.*;

int NUM_ITERATIONS = 10;

GrayScott gs;
ToneMap toneMap;

void setup() {
  size(256,256);
  gs=new GrayScott(width,height,false);
  gs.setCoefficients(0.021,0.076,0.12,0.06);
  // create a color gradient for 256 values
  ColorGradient grad=new ColorGradient();
  // NamedColors are preset colors, but any TColor can be added
  // see javadocs for list of names:
  // http://toxiclibs.org/docs/colorutils/toxi/color/NamedColor.html
  grad.addColorAt(0,NamedColor.BLACK);
  grad.addColorAt(128,NamedColor.RED);
  grad.addColorAt(192,NamedColor.YELLOW);
  grad.addColorAt(255,NamedColor.WHITE);
  // this gradient is used to map simulation values to colors
  // the first 2 parameters define the min/max values of the
  // input range (Gray-Scott produces values in the interval of 0.0 - 0.5)
  // setting the max = 0.33 increases the contrast
  toneMap=new ToneMap(0,0.33,grad);
}

void draw() {
  if (mousePressed) {
    gs.setRect(mouseX, mouseY,20,20);
  }
  loadPixels();
  // update the simulation a few time steps
  for(int i=0; i<NUM_ITERATIONS; i++) {
    gs.update(1);
  }
  // read out the V result array
  // and use tone map to render colours
  for(int i=0; i<gs.v.length; i++) {
    pixels[i]=toneMap.getARGBToneFor(gs.v[i]);
  }
  updatePixels();
}

void keyPressed() {
  gs.reset();
}