Search in sources :

Example 6 with Vector3D

use of mcib3d.geom.Vector3D in project mcib3d-core by mcib3d.

the class DeformableMesh method getGlobalScaling.

public double getGlobalScaling() {
    Point3D center = this.getCenterAsVector();
    double ratio = 0;
    double r;
    int cpt = 0;
    for (int i = 0; i < forces.size(); i++) {
        Point3D P0 = new Point3D(this.getUniqueVertex(i));
        Point3D P1 = new Point3D(this.getUniqueVertex(i));
        Vector3D force = forces.get(i);
        // IJ.log("Force " + i + " " + force);
        if (force != null) {
            P1.translate(force);
            Vector3D V0 = new Vector3D(center, P0);
            Vector3D V1 = new Vector3D(center, P1);
            r = V1.getLength() / V0.getLength();
            ratio += r;
            cpt++;
        // IJ.log("scale " + i + " " + r + " " + V0 + " " + V1);
        }
    }
    if (cpt > 1) {
        ratio /= (double) cpt;
    } else {
        ratio = 1;
    }
    return ratio;
}
Also used : Vector3D(mcib3d.geom.Vector3D) Point3D(mcib3d.geom.Point3D)

Example 7 with Vector3D

use of mcib3d.geom.Vector3D in project mcib3d-core by mcib3d.

the class DeformableMesh method getRotationAxis.

public Vector3D getRotationAxis() {
    Vector3D moy = new Vector3D();
    Vector3D elong = this.getMainAxis().getNormalizedVector();
    double cpt = 0;
    for (int i = 0; i < getNbUniqueVertices(); i++) {
        Vector3D displ = forces.get(i);
        if (displ != null) {
            // double dot=displN.dotProduct(elong);
            if (displ.getX() < 0) {
                displ = displ.multiply(-1);
            }
            moy.addMe(displ.crossProduct(elong));
            cpt++;
        }
    }
    if (cpt > 0) {
        moy = moy.multiply(1.0 / cpt);
    }
    return moy.getNormalizedVector();
}
Also used : Vector3D(mcib3d.geom.Vector3D)

Example 8 with Vector3D

use of mcib3d.geom.Vector3D in project mcib3d-core by mcib3d.

the class DeformableMesh method getTranslation.

public Vector3D getTranslation() {
    Vector3D moy = new Vector3D();
    int cpt = 0;
    for (int i = 0; i < getNbUniqueVertices(); i++) {
        Vector3D displ = forces.get(i);
        if (displ != null) {
            moy.addMe(displ);
            cpt++;
        }
    }
    if (cpt > 0) {
        moy = moy.multiply(1.0 / (double) (cpt));
    }
    return moy;
}
Also used : Vector3D(mcib3d.geom.Vector3D)

Example 9 with Vector3D

use of mcib3d.geom.Vector3D in project mcib3d-core by mcib3d.

the class SymmetryFilter method computeBins.

private void computeBins() {
    ImageHandler img = edges[0];
    bin1 = new ImageFloat("bin1", img.sizeX, img.sizeY, img.sizeZ);
    bin2 = new ImageFloat("bin2", img.sizeX, img.sizeY, img.sizeZ);
    for (int z = 0; z < img.sizeZ; z++) {
        IJ.showStatus("Symmetry " + z + "/" + img.sizeZ);
        for (int x = 0; x < img.sizeX; x++) {
            for (int y = 0; y < img.sizeY; y++) {
                double ex = edges[0].getPixel(x, y, z);
                double ey = edges[1].getPixel(x, y, z);
                double ez = edges[2].getPixel(x, y, z);
                double ee = Math.sqrt(ex * ex + ey * ey + ez * ez);
                // bin
                Vector3D grad = new Vector3D(ex, ey, ez);
                grad.normalize();
                if (grad.getLength() == 0) {
                    continue;
                }
                Point3D pos = new Vector3D(x, y, z);
                for (int d = 0; d < radius; d++) {
                    pos.translate(grad);
                    if ((d > 0) && (img.contains(pos.getRoundX(), pos.getRoundY(), pos.getRoundZ()))) {
                        bin1.setPixelIncrement(pos, 1);
                        if (improved) {
                            bin2.setPixelIncrement(pos, (float) (d * ee));
                        } else {
                            bin2.setPixelIncrement(pos, (float) (ee));
                        }
                    }
                }
            }
        }
    }
}
Also used : ImageHandler(mcib3d.image3d.ImageHandler) Vector3D(mcib3d.geom.Vector3D) Point3D(mcib3d.geom.Point3D) ImageFloat(mcib3d.image3d.ImageFloat)

Example 10 with Vector3D

use of mcib3d.geom.Vector3D in project mcib3d-core by mcib3d.

the class FHTImage3D method fill3D.

/**
 * Fill 3D Reconstruction using central section theorem
 *
 * @param proj Array of FFT of the projections
 * @param W    Orientations of the images
 */
public void fill3D(FHTImage3D[] proj, Vector3D[] W) {
    if (!frequencyDomain) {
        return;
    }
    Vector3D Y = new Vector3D(0, 1, 0);
    int sizexy = sizex * sizey;
    double ox = sizex / 2;
    double oy = sizey / 2;
    double oz = sizez / 2;
    double x0;
    double y0;
    double z0;
    double xc0;
    double yc0;
    double zc0;
    int xx0;
    int yy0;
    int nbproj = proj.length;
    double[] a11 = new double[nbproj];
    double[] a12 = new double[nbproj];
    double[] a21 = new double[nbproj];
    double[] a22 = new double[nbproj];
    double[] a31 = new double[nbproj];
    double[] a32 = new double[nbproj];
    double[] at11 = new double[nbproj];
    double[] at12 = new double[nbproj];
    double[] at13 = new double[nbproj];
    double[] at21 = new double[nbproj];
    double[] at22 = new double[nbproj];
    double[] at23 = new double[nbproj];
    Vector3D U;
    Vector3D V;
    double RP0;
    double XX0;
    double XX1;
    double XX2;
    double RP1;
    for (int p = 0; p < nbproj; p++) {
        U = Y.crossProduct(W[p]);
        V = W[p].crossProduct(U);
        U.normalize();
        V.normalize();
        double ux = U.getX();
        double uy = U.getY();
        double uz = U.getZ();
        double vx = V.getX();
        double vy = V.getY();
        double vz = V.getZ();
        double wx = W[p].getX();
        double wy = W[p].getY();
        double wz = W[p].getZ();
        // A = new Matrix(3, 3);
        // A.set(0, 0, ux);
        // A.set(0, 1, vx);
        // A.set(0, 2, wx);
        // A.set(1, 0, uy);
        // A.set(1, 1, vy);
        // A.set(1, 2, wy);
        // A.set(2, 0, uz);
        // A.set(2, 1, vz);
        // A.set(2, 2, wz);
        a11[p] = ux;
        a12[p] = vx;
        a21[p] = uy;
        a22[p] = vy;
        a31[p] = uz;
        a32[p] = vz;
        // AT = A.inverse();
        // at11[p] = AT.get(0, 0);
        // at12[p] = AT.get(0, 1);
        // at13[p] = AT.get(0, 2);
        // at21[p] = AT.get(1, 0);
        // at22[p] = AT.get(1, 1);
        // at23[p] = AT.get(1, 2);
        at11[p] = ux;
        at12[p] = uy;
        at13[p] = uz;
        at21[p] = vx;
        at22[p] = vy;
        at23[p] = vz;
    }
    double dist;
    double distmin;
    int projmin;
    Chrono time = new Chrono(sizex);
    time.start();
    for (int x = 0; x < sizex; x++) {
        for (int y = 0; y < sizey; y++) {
            for (int z = 0; z < sizez; z++) {
                x0 = x - ox;
                y0 = y - oy;
                z0 = z - oz;
                distmin = 1.0;
                projmin = -1;
                for (int p = 0; p < nbproj; p++) {
                    RP0 = at11[p] * x0 + at12[p] * y0 + at13[p] * z0;
                    RP1 = at21[p] * x0 + at22[p] * y0 + at23[p] * z0;
                    xx0 = (int) Math.round(RP0);
                    yy0 = (int) Math.round(RP1);
                    XX0 = a11[p] * xx0 + a12[p] * yy0;
                    XX1 = a21[p] * xx0 + a22[p] * yy0;
                    XX2 = a31[p] * xx0 + a32[p] * yy0;
                    xc0 = XX0 + ox;
                    yc0 = XX1 + oy;
                    zc0 = XX2 + oz;
                    dist = (x - xc0) * (x - xc0) + (y - yc0) * (y - yc0) + (z - zc0) * (z - zc0);
                    if (dist <= distmin) {
                        distmin = dist;
                        projmin = p;
                    }
                }
                // index = 0;
                if (projmin != -1) {
                    for (int k = -1; k <= 1; k++) {
                        for (int l = -1; l <= 1; l++) {
                            RP0 = at11[projmin] * x0 + at12[projmin] * y0 + at13[projmin] * z0;
                            RP1 = at21[projmin] * x0 + at22[projmin] * y0 + at23[projmin] * z0;
                            xx0 = (int) Math.round(RP0 + k);
                            yy0 = (int) Math.round(RP1 + l);
                            XX0 = a11[projmin] * xx0 + a12[projmin] * yy0;
                            XX1 = a21[projmin] * xx0 + a22[projmin] * yy0;
                            XX2 = a31[projmin] * xx0 + a32[projmin] * yy0;
                            xc0 = XX0 + ox;
                            yc0 = XX1 + oy;
                            zc0 = XX2 + oz;
                            // ll[index] = W[projmin].intersection_unit_cube(x, y, z, xc0, yc0, zc0);
                            // real = proj[projmin].getPixelReal((int) (xx0 + ox), (int) (yy0 + oy), 0);
                            // imag = proj[projmin].getPixelImag((int) (xx0 + ox), (int) (yy0 + oy), 0);
                            // this.putPixel(x, y, z, getPixelReal(x, y, z) + (float) real, getPixelImag(x, y, z) + (float) imag);
                            xx0 += ox;
                            yy0 += oy;
                            // }
                            if (xx0 > 0 && xx0 < sizex && yy0 > 0 && yy0 < sizey) {
                                int fhtcoord = (int) (xx0 + yy0 * sizex);
                                putPixel(x, y, z, getPixel(x, y, z) + proj[projmin].getPixel(xx0, yy0, 0));
                            }
                        // index++;
                        }
                    }
                }
            }
        }
        time.stop();
        System.out.print("\r                                                                 \r" + (100 * (x + 1) / sizex) + "% \t" + time.delayString() + "\t (" + time.remainString(x + 1) + ")           ");
    }
/*
         *  /////////////////////////////////////////////////////
         *  int nbproj = proj.length;
         *  Matrix[] M = new Matrix[nbproj];
         *  Matrix MM;
         *  Vector3D U;
         *  Vector3D V;
         *  /Vector3D X = new Vector3D(0, 0, 1);
         *  Vector3D Y = new Vector3D(0, 1, 0);
         *  Matrix B;
         *  Matrix res;
         *  double[] b = new double[6];
         *  double dist;
         *  int projmin;
         *  double amin;
         *  double bmin;
         *  double distmin;
         *  float real;
         *  float imag;
         *  float coeff;
         *  double alpha;
         *  double beta;
         *  double ux[] = new double[nbproj];
         *  double uy[] = new double[nbproj];
         *  double uz[] = new double[nbproj];
         *  double vx[] = new double[nbproj];
         *  double vy[] = new double[nbproj];
         *  double vz[] = new double[nbproj];
         *  double ox = 0;
         *  double oy = 0;
         *  double oz = 0;
         *  double wx[] = new double[nbproj];
         *  double wy[] = new double[nbproj];
         *  double wz[] = new double[nbproj];
         *  double ww[] = new double[nbproj];
         *  for (int p = 0; p < nbproj; p++) {
         *  U = Y.crossProduct(W[p]);
         *  V = U.crossProduct(W[p]);
         *  U.normalize();
         *  V.normalize();
         *  ux[p] = U.getX();
         *  uy[p] = U.getY();
         *  uz[p] = U.getZ();
         *  vx[p] = V.getX();
         *  vy[p] = V.getY();
         *  vz[p] = V.getZ();
         *  ox = getSizex() / 2;
         *  oy = getSizey() / 2;
         *  oz = getSizez() / 2;
         *  wx[p] = W[p].getX();
         *  wy[p] = W[p].getY();
         *  wz[p] = W[p].getZ();
         *  ww[p] = wx[p] * wx[p] + wy[p] * wy[p] + wz[p] * wz[p];
         *  }
         *  int sx = getSizex();
         *  int sy = getSizey();
         *  int sz = getSizez();
         *  double px;
         *  double py;
         *  double pz;
         *  double lambda;
         *  for (int x = 0; x < sx; x++) {
         *  for (int y = 0; y < sy; y++) {
         *  for (int z = 0; z < sz; z++) {
         *  / parcours des projections
         *  distmin = 1.0;
         *  projmin = -1;
         *  coeff = 0;
         *  amin = 0;
         *  bmin = 0;
         *  real = 0;
         *  imag = 0;
         *  double xxmin = 0;
         *  double yymin = 0;
         *  for (int p = 0; p < nbproj; p++) {
         *  / resolution directe
         *  lambda = (wx[p] * (ox - x) + wy[p] * (oy - y) + wz[p] * (oz - z)) / ww[p];
         *  px = x + lambda * wx[p];
         *  py = y + lambda * wy[p];
         *  pz = z + lambda * wz[p];
         *  alpha = (px - ox) * ux[p] + (py - oy) * uy[p] + (pz - oz) * uz[p];
         *  beta = (px - ox) * vx[p] + (py - oy) * vy[p] + (pz - oz) * vz[p];
         *  double xx = (alpha + sx / 2);
         *  double yy = (sy / 2 - beta);
         *  dist = Math.abs(lambda);
         *  if (dist < distmin) {
         *  distmin = dist;
         *  projmin = p;
         *  xxmin = xx;
         *  yymin = yy;
         *  }
         *  }
         *  if (projmin != -1) {
         *  real = proj[projmin].getInterpolatedPixelReal((float) xxmin, (float) yymin, 0);
         *  imag = proj[projmin].getInterpolatedPixelImag((float) xxmin, (float) yymin, 0);
         *  this.putPixel(x, y, z, real, imag);
         *  }
         *  }
         *  }
         *  }
         */
}
Also used : Chrono(mcib3d.utils.Chrono) Vector3D(mcib3d.geom.Vector3D)

Aggregations

Vector3D (mcib3d.geom.Vector3D)18 Point3D (mcib3d.geom.Point3D)4 EigenvalueDecomposition (mcib3d.Jama.EigenvalueDecomposition)1 Matrix (mcib3d.Jama.Matrix)1 Object3D (mcib3d.geom.Object3D)1 Objects3DPopulation (mcib3d.geom.Objects3DPopulation)1 ImageFloat (mcib3d.image3d.ImageFloat)1 ImageHandler (mcib3d.image3d.ImageHandler)1 ArrayUtil (mcib3d.utils.ArrayUtil)1 Chrono (mcib3d.utils.Chrono)1 Point3f (org.scijava.vecmath.Point3f)1