Examples with PoissonDistribution org.apache.commons.math3.distribution.PoissonDistribution used on opensource projects

Example 16 with PoissonDistribution

use of org.apache.commons.math3.distribution.PoissonDistribution in project vcell by virtualcell.

the class BrownianDynamics2DSolver method sampleImage.

public UShortImage sampleImage(int numX, int numY, Origin origin, Extent extent, double diffusion, double integrationTime, int numSteps, double psfVar, double detectorGainAndQuantumYeild, boolean bConvolve, boolean bNoise) throws ImageException {
    double deltaT = integrationTime / (numSteps);
    // 
    // sample the x and y coordinates for the structured grid (pixel coordinates). .. coordinates will be used lots of times (let's precompute).
    // 
    double[] imageX = new double[numX];
    double[] imageY = new double[numY];
    for (int i = 0; i < numX; i++) {
        imageX[i] = origin.getX() + (0.5 + i) * extent.getX() / numX;
    }
    for (int j = 0; j < numY; j++) {
        imageY[j] = origin.getY() + (0.5 + j) * extent.getY() / numY;
    }
    // apply point spread function to gather poisson rates from all particles to all pixels.
    // then we roll the dice to sample from poission later.
    double[] particleImage = new double[numX * numY];
    Arrays.fill(particleImage, 0);
    for (int step = 0; step < numSteps; step++) {
        if (bConvolve) {
            // 
            // convolve all particles with normalized psf (each particle integrates to one over neighboring pixels) for each pixel and accumulate in lambda array
            // (lambda is not reset so we can repeat and integrate over time)
            // 
            System.out.println("convolving particles, step=" + step + ", time=" + currTime);
            addConvolved(numX, numY, psfVar, imageX, imageY, particleImage);
        } else {
            // 
            // count particles in each pixel
            // 
            System.out.println("binning particles, step=" + step + ", time=" + currTime);
            addBinned(numX, numY, origin.getX(), origin.getY(), extent.getX(), extent.getY(), particleImage);
        }
        // 
        // let diffusion happen during sampling
        // 
        System.out.println("simulating particles, step=" + step + ", time=" + currTime);
        step(diffusion, deltaT);
    }
    // 
    // generate sample image by generating fluorescence intensity
    // need to rescale to consider detector gain and expected number of photons from each particle. (ignore detector noise)
    // 
    float[] pixels = new float[particleImage.length];
    double totalGain = detectorGainAndQuantumYeild * integrationTime / numSteps;
    for (int i = 0; i < pixels.length; i++) {
        float rate = (float) (particleImage[i] * totalGain);
        if (rate > 1e-5) {
            if (bNoise) {
                pixels[i] = new PoissonDistribution(imagingRng, rate, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS).sample();
            } else {
                pixels[i] = rate;
            }
        }
    }
    // 
    // create unsigned short sample image
    // 
    short[] ushortPixels = new short[pixels.length];
    for (int i = 0; i < ushortPixels.length; i++) {
        if (pixels[i] > 32767) {
            throw new RuntimeException("pixel overflow of signed short, value = " + pixels[i]);
        }
        if (pixels[i] < 0) {
            throw new RuntimeException("unexpected negative pixel, value = " + pixels[i]);
        }
        ushortPixels[i] = (short) (0xffff & (int) pixels[i]);
    }
    UShortImage sampledImage = new UShortImage(ushortPixels, origin, extent, numX, numY, 1);
    return sampledImage;
}

Example 17 with PoissonDistribution

use of org.apache.commons.math3.distribution.PoissonDistribution in project vcell by virtualcell.

the class BrownianDynamics2DSolver method bleachGuassian0.

private long bleachGuassian0(double muX, double muY, double psfVar, double bleachFactor_K, double deltaT) {
    // 
    // distribute each particle using a normalized Gaussian point spread function with variance psfVar.
    // we add the fluorescence to the image (can be called multiple times to accumulate fluorescence).
    // 
    long bleachCount = 0;
    double DISTANCE_6_SIGMA = Math.sqrt(psfVar) * 6;
    for (int p = 0; p < numParticles; p++) {
        if (bFluorescent[p]) {
            double pX = particleX[p] - muX;
            double pY = particleY[p] - muY;
            double radius2 = pX * pX + pY * pY;
            if (radius2 < DISTANCE_6_SIGMA) {
                double particleFluorRate = bleachFactor_K * deltaT * FastMath.exp(-radius2 / psfVar);
                PoissonDistribution poisson = new PoissonDistribution(simulationRng, particleFluorRate, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
                if (poisson.sample() >= 1) {
                    // if bleaching event happened
                    bFluorescent[p] = false;
                    bleachCount++;
                }
            }
        }
    }
    return bleachCount;
}