Examples with Lightsheet net.preibisch.simulation.raytracing.Lightsheet used on opensource projects

Example 1 with Lightsheet

use of net.preibisch.simulation.raytracing.Lightsheet in project multiview-simulation by PreibischLab.

the class SimulateMultiViewAberrations method refract3d.

public static VolumeInjection refract3d(final RandomAccessibleInterval<FloatType> imgIn, final RandomAccessibleInterval<FloatType> imgRi, final boolean illum, final int z, final double lsMiddle, final double lsEdge, final double ri) {
    // the refracted image
    final Img<FloatType> img = new ArrayImgFactory<FloatType>(new FloatType()).create(imgIn);
    final Img<FloatType> weight = new ArrayImgFactory<FloatType>(new FloatType()).create(imgIn);
    // for drawing individual rays ...
    // final RandomAccess< FloatType > rOut = img.randomAccess();
    // final long[] intPos = new long[ 3 ];
    final RealRandomAccess<FloatType> rrIm = Views.interpolate(Views.extendMirrorSingle(imgIn), new NLinearInterpolatorFactory<>()).realRandomAccess();
    final RealRandomAccess<FloatType> rrRi = Views.interpolate(Views.extendMirrorSingle(imgRi), new NLinearInterpolatorFactory<>()).realRandomAccess();
    // row, column
    final double[][] matrix = new double[3][3];
    final double[] eigenVector = new double[3];
    final double[] sigma = new double[] { 0.5, 0.5, 0.5 };
    final VolumeInjection inject = new VolumeInjection(img, weight, sigma);
    // the incoming direction of the light
    final double[] rayVector = new double[3];
    // the outgoing, refracted direction of the light
    final double[] refractedRay = new double[3];
    // the current position of the ray
    final double[] rayPosition = new double[3];
    // air (intensity == 0)
    final double nA = 1.00;
    // water (intensity == 1)
    final double nB = ri;
    System.out.println("middle: " + lsMiddle);
    System.out.println("lsEdge: " + lsEdge);
    final Lightsheet ls = new Lightsheet(img.dimension(0) / 2.0, lsMiddle, img.dimension(0), lsEdge);
    System.out.println("starting...");
    long time = System.currentTimeMillis();
    // for each point on the xz plane at y=0
    final int numRays = 200000;
    final long maxMoves = img.dimension(2);
    final Random rnd = new Random(2423);
    for (int i = 0; i < numRays; ++i) {
        // x;//c.getIntPosition( 0 );
        rayPosition[0] = rnd.nextDouble() * imgIn.max(0);
        rayPosition[1] = illum ? (int) imgIn.dimension(1) - 1 : 0;
        final double lightsheetthickness = ls.predict(rayPosition[0]);
        // c.getIntPosition( 1 );
        rayPosition[2] = z + (rnd.nextDouble() * lightsheetthickness) - lightsheetthickness / 2.0;
        rayVector[0] = (rnd.nextDouble() - 0.5) / 5;
        rayVector[1] = illum ? -1 : 1;
        rayVector[2] = 0;
        Raytrace.norm(rayVector);
        int moves = 0;
        while (inside(rayPosition, imgIn) && moves < maxMoves) {
            ++moves;
            rrIm.setPosition(rayPosition);
            rrRi.setPosition(rayPosition);
            final float valueIm = rrIm.get().get();
            // normal vector of refraction plane (still maybe needs to be inverted to point towards the incoming signal)
            Hessian.computeHessianMatrix3D(rrRi, matrix);
            double ev = Hessian.computeLargestEigenVectorAndValue3d(matrix, eigenVector);
            if (Math.abs(ev) > 0.01) {
                // compute refractive index change
                rrRi.setPosition(rayPosition);
                rrRi.move(-rayVector[0], 0);
                rrRi.move(-rayVector[1], 1);
                rrRi.move(-rayVector[2], 2);
                // intensity at the origin of the ray
                final double i0 = rrRi.get().get();
                rrRi.move(2 * rayVector[0], 0);
                rrRi.move(2 * rayVector[1], 1);
                rrRi.move(2 * rayVector[2], 2);
                // intensity at the projected location of the ray
                final double i1 = rrRi.get().get();
                final double n0 = (nB - nA) * i0 + nA;
                final double n1 = (nB - nA) * i1 + nA;
                final double thetaI = Raytrace.incidentAngle(rayVector, eigenVector);
                final double thetaT = Raytrace.refract(rayVector, eigenVector, n0, n1, thetaI, refractedRay);
                // total reflection
                if (Double.isNaN(thetaT)) {
                    refractedRay[0] = rayVector[0];
                    refractedRay[1] = rayVector[1];
                    refractedRay[2] = rayVector[2];
                // continue; // that's an expensive getting stuck
                // Raytrace.reflect( rayVector, eigenVector, refractedRay );
                }
                Raytrace.norm(refractedRay);
                // update the ray vector
                rayVector[0] = refractedRay[0];
                rayVector[1] = refractedRay[1];
                rayVector[2] = refractedRay[2];
            }
            // place a gaussian sphere
            inject.addNormalizedGaussian(valueIm, rayPosition);
            /*
				// for drawing individual rays ...
				intPos[ 0 ] = Math.round( rayPosition[ 0 ] );
				intPos[ 1 ] = Math.round( rayPosition[ 1 ] );
				intPos[ 2 ] = Math.round( rayPosition[ 2 ] );
				rOut.setPosition( intPos );
				rOut.get().set( value );
				*/
            rayPosition[0] += rayVector[0];
            rayPosition[1] += rayVector[1];
            rayPosition[2] += rayVector[2];
        }
    }
    System.out.println(" ... " + (System.currentTimeMillis() - time));
    return inject;
}