Examples with TaubinSmoothing cbit.vcell.geometry.surface.TaubinSmoothing used on opensource projects

Example 1 with TaubinSmoothing

use of cbit.vcell.geometry.surface.TaubinSmoothing in project vcell by virtualcell.

the class PDEDataViewer method createDataValueSurfaceViewer.

/**
 * Insert the method's description here.
 * Creation date: (9/25/2005 1:53:00 PM)
 */
private DataValueSurfaceViewer createDataValueSurfaceViewer(ClientTaskStatusSupport clientTaskStatusSupport) throws ImageException, UserCancelException {
    // try{
    // if(fieldDataValueSurfaceViewer == null){
    // Surfaces
    CartesianMesh cartesianMesh = getPdeDataContext().getCartesianMesh();
    if (cartesianMesh.getMembraneElements() == null || cartesianMesh.getMembraneElements().length == 0 || cartesianMesh.isChomboMesh()) {
        // return fieldDataValueSurfaceViewer;
        return getDataValueSurfaceViewer();
    }
    meshRegionSurfaces = new MeshDisplayAdapter(cartesianMesh).generateMeshRegionSurfaces(clientTaskStatusSupport);
    SurfaceCollection surfaceCollection = meshRegionSurfaces.getSurfaceCollection();
    // SurfaceNames
    final String[] surfaceNames = new String[meshRegionSurfaces.getSurfaceCollection().getSurfaceCount()];
    for (int i = 0; i < meshRegionSurfaces.getSurfaceCollection().getSurfaceCount(); i++) {
        // Get the first element, any will do, all have same inside/outside volumeIndex
        MembraneElement me = cartesianMesh.getMembraneElements()[meshRegionSurfaces.getMembraneIndexForPolygon(i, 0)];
        if (getSimulationModelInfo() != null) {
            surfaceNames[i] = getSimulationModelInfo().getMembraneName(cartesianMesh.getSubVolumeFromVolumeIndex(me.getInsideVolumeIndex()), cartesianMesh.getSubVolumeFromVolumeIndex(me.getOutsideVolumeIndex()), false);
        } else {
            surfaceNames[i] = i + "";
        }
    }
    // SurfaceAreas
    final Double[] surfaceAreas = new Double[meshRegionSurfaces.getSurfaceCollection().getSurfaceCount()];
    for (int i = 0; i < meshRegionSurfaces.getSurfaceCollection().getSurfaceCount(); i++) {
        surfaceAreas[i] = new Double(cartesianMesh.getRegionMembraneSurfaceAreaFromMembraneIndex(meshRegionSurfaces.getMembraneIndexForPolygon(i, 0)));
    }
    // DataValueSurfaceViewer fieldDataValueSurfaceViewer0 = new DataValueSurfaceViewer();
    TaubinSmoothing taubinSmoothing = new TaubinSmoothingWrong();
    TaubinSmoothingSpecification taubinSpec = TaubinSmoothingSpecification.getInstance(.3);
    taubinSmoothing.smooth(surfaceCollection, taubinSpec);
    getDataValueSurfaceViewer().init(meshRegionSurfaces.getSurfaceCollection(), cartesianMesh.getOrigin(), cartesianMesh.getExtent(), surfaceNames, surfaceAreas, cartesianMesh.getGeometryDimension());
    return getDataValueSurfaceViewer();
// }
// }catch(UserCancelException e){
// throw e;
// }catch(Exception e){
// PopupGenerator.showErrorDialog(PDEDataViewer.this, e.getMessage(), e);
// }
// return fieldDataValueSurfaceViewer;
}

Example 2 with TaubinSmoothing

use of cbit.vcell.geometry.surface.TaubinSmoothing in project vcell by virtualcell.

the class RegionImage method calculateRegions_New.

// 
// Calculate regions using single pass algorithm.  Creates information
// used to generate surfaces as well
// 
private void calculateRegions_New(VCImage vcImage, int dimension, Extent extent, Origin origin, ClientTaskStatusSupport clientTaskStatusSupport) throws ImageException {
    // Find linked pixel values in x,y,z and surface elements locations
    final int EMPTY_REGION = 0;
    final int FIRST_REGION = 1;
    final int ARRAY_SIZE_INCREMENT = 1000;
    byte[] imagePixels = vcImage.getPixels();
    int[] regionPixels = new int[imagePixels.length];
    BitSet xSurfElements = new BitSet(imagePixels.length);
    BitSet ySurfElements = new BitSet(imagePixels.length);
    BitSet zSurfElements = new BitSet(imagePixels.length);
    Vector<Integer>[] regionLinkArr = new Vector[ARRAY_SIZE_INCREMENT];
    int[] regionSizeArr = new int[ARRAY_SIZE_INCREMENT];
    Vector<Byte> regionImagePixelV = new Vector<Byte>();
    // 0 not used
    regionImagePixelV.add((byte) 0);
    byte currentImagePixel;
    int currentRegion;
    int masterIndex = 0;
    int nextAvailableRegion = FIRST_REGION;
    for (int zIndex = 0; zIndex < vcImage.getNumZ(); zIndex++) {
        if (clientTaskStatusSupport != null) {
            if (clientTaskStatusSupport.isInterrupted()) {
                return;
            }
        }
        int zForwardIndex = ((zIndex + 1) < vcImage.getNumZ() ? masterIndex + (vcImage.getNumX() * vcImage.getNumY()) : -1);
        for (int yIndex = 0; yIndex < vcImage.getNumY(); yIndex++) {
            int yForwardIndex = ((yIndex + 1) < vcImage.getNumY() ? masterIndex + vcImage.getNumX() : -1);
            currentImagePixel = imagePixels[masterIndex];
            if (regionPixels[masterIndex] != EMPTY_REGION) {
                currentRegion = regionPixels[masterIndex];
            } else {
                currentRegion = nextAvailableRegion;
                if (currentRegion >= regionLinkArr.length) {
                    // make more room for arrays
                    Vector<Integer>[] temp = new Vector[regionLinkArr.length + ARRAY_SIZE_INCREMENT];
                    System.arraycopy(regionLinkArr, 0, temp, 0, regionLinkArr.length);
                    regionLinkArr = temp;
                    int[] temp2 = new int[regionSizeArr.length + ARRAY_SIZE_INCREMENT];
                    System.arraycopy(regionSizeArr, 0, temp2, 0, regionSizeArr.length);
                    regionSizeArr = temp2;
                }
                regionPixels[masterIndex] = currentRegion;
                if (regionImagePixelV.size() != currentRegion) {
                    throw new ImageException("Mismatch between region and pixel buffer");
                }
                regionImagePixelV.add(currentImagePixel);
                nextAvailableRegion += 1;
            }
            for (int xIndex = 0; xIndex < vcImage.getNumX(); xIndex++) {
                if (imagePixels[masterIndex] == currentImagePixel) {
                    if (regionPixels[masterIndex] != EMPTY_REGION) {
                        if (currentRegion != regionPixels[masterIndex]) {
                            createLink(regionLinkArr, currentRegion, regionPixels, masterIndex);
                        }
                    } else {
                        regionPixels[masterIndex] = currentRegion;
                    }
                } else {
                    xSurfElements.set(masterIndex - 1);
                    currentImagePixel = imagePixels[masterIndex];
                    if (regionPixels[masterIndex] != EMPTY_REGION) {
                        currentRegion = regionPixels[masterIndex];
                    } else {
                        currentRegion = nextAvailableRegion;
                        if (currentRegion >= regionLinkArr.length) {
                            // make more room for arrays
                            Vector<Integer>[] temp = new Vector[regionLinkArr.length + ARRAY_SIZE_INCREMENT];
                            System.arraycopy(regionLinkArr, 0, temp, 0, regionLinkArr.length);
                            regionLinkArr = temp;
                            int[] temp2 = new int[regionSizeArr.length + ARRAY_SIZE_INCREMENT];
                            System.arraycopy(regionSizeArr, 0, temp2, 0, regionSizeArr.length);
                            regionSizeArr = temp2;
                        }
                        regionPixels[masterIndex] = currentRegion;
                        if (regionImagePixelV.size() != currentRegion) {
                            throw new ImageException("Mismatch between region and pixel buffer");
                        }
                        regionImagePixelV.add(currentImagePixel);
                        nextAvailableRegion += 1;
                    }
                }
                regionSizeArr[currentRegion] += 1;
                // }
                if (yForwardIndex != -1) {
                    if (imagePixels[yForwardIndex] == currentImagePixel) {
                        if (regionPixels[yForwardIndex] == EMPTY_REGION) {
                            regionPixels[yForwardIndex] = currentRegion;
                        } else if (currentRegion != regionPixels[yForwardIndex]) {
                            createLink(regionLinkArr, currentRegion, regionPixels, yForwardIndex);
                        }
                    } else {
                        ySurfElements.set(masterIndex);
                    }
                    yForwardIndex += 1;
                }
                if (zForwardIndex != -1) {
                    if (imagePixels[zForwardIndex] == currentImagePixel) {
                        if (regionPixels[zForwardIndex] == EMPTY_REGION) {
                            regionPixels[zForwardIndex] = currentRegion;
                        } else if (currentRegion != regionPixels[zForwardIndex]) {
                            createLink(regionLinkArr, currentRegion, regionPixels, zForwardIndex);
                        }
                    } else {
                        zSurfElements.set(masterIndex);
                    }
                    zForwardIndex += 1;
                }
                masterIndex += 1;
            }
        }
    }
    // System.out.println(xSurfElements.cardinality()+" "+ySurfElements.cardinality()+" "+zSurfElements.cardinality());
    // System.out.println("----------link time "+((System.currentTimeMillis()-startTime)/1000.0));
    // startTime = System.currentTimeMillis();
    // Distribute links
    Vector<Integer>[] collector = new Vector[nextAvailableRegion];
    for (int i = 1; i < nextAvailableRegion; /*regionSizeArr.length*/
    i++) {
        // 0 not used
        if (regionSizeArr[i] != 0) {
            collector[i] = (Vector) (regionLinkArr[i] != null ? regionLinkArr[i].clone() : null);
        }
    }
    for (int i = 1; i < nextAvailableRegion; /*regionSizeArr.length*/
    i++) {
        // 0 not used
        if (clientTaskStatusSupport != null && clientTaskStatusSupport.isInterrupted()) {
            return;
        }
        // System.out.print("region="+i+" size="+regionSizeArr[i]);
        Vector<Integer> intV = regionLinkArr[i];
        for (int j = 0; intV != null && j < intV.size(); j++) {
            // System.out.print((j==0?" - ":" ")+intV.elementAt(j));
            // Collect
            Vector<Integer> collectIntV = collector[intV.elementAt(j)];
            if (collectIntV == null) {
                collectIntV = new Vector<Integer>();
                collector[intV.elementAt(j)] = collectIntV;
                if (regionLinkArr[intV.elementAt(j)] != null) {
                    collectIntV.addAll(regionLinkArr[intV.elementAt(j)]);
                }
            }
            if (!collectIntV.contains(i)) {
                collectIntV.add(i);
            }
        }
    // System.out.println();
    }
    // System.out.println("----------distribute link time "+((System.currentTimeMillis()-startTime)/1000.0));
    // startTime = System.currentTimeMillis();
    // for (int i = 1; i < collector.length; i++) {// 0 not used
    // if(regionSizeArr[i] != 0){
    // System.out.print("Collect region="+i+" size="+regionSizeArr[i]);
    // Vector<Integer> intV = collector[i];
    // for (int j = 0; intV!= null && j < intV.size();j++) {
    // System.out.print((j==0?" - ":" ")+intV.elementAt(j));
    // }
    // System.out.println();
    // }
    // }
    // Gather links into distinct regions
    // Map link-regions(implicit) to distinct-regions
    int[] linkRegionMap = new int[collector.length];
    int totalSize = 0;
    Vector<Vector<Integer>> regionsV = new Vector<Vector<Integer>>();
    Vector<Integer> regionsSizeV = new Vector<Integer>();
    BitSet checkFlagBS = new BitSet(collector.length);
    for (int i = 1; i < collector.length; i++) {
        // 0 not used
        if (clientTaskStatusSupport != null && clientTaskStatusSupport.isInterrupted()) {
            return;
        }
        if (checkFlagBS.get(i)) {
            continue;
        }
        checkFlagBS.set(i);
        Vector<Integer> holderV = new Vector<Integer>();
        holderV.add(i);
        if (collector[i] != null) {
            holderV.addAll(collector[i]);
        }
        int checkIndex = 0;
        // 
        boolean[] holderVContainsFlag = new boolean[collector.length];
        for (int j = 0; j < holderV.size(); j++) {
            holderVContainsFlag[holderV.elementAt(j)] = true;
        }
        // 
        while (true) {
            if (collector[holderV.elementAt(checkIndex)] != null) {
                Vector<Integer> newLinksV = collector[holderV.elementAt(checkIndex)];
                for (int j = 0; j < newLinksV.size(); j++) {
                    if (!checkFlagBS.get(newLinksV.elementAt(j)) && !holderVContainsFlag[newLinksV.elementAt(j)]) // !holderV.contains(newLinksV.elementAt(j))
                    {
                        holderV.add(newLinksV.elementAt(j));
                        holderVContainsFlag[newLinksV.elementAt(j)] = true;
                    }
                }
            }
            checkFlagBS.set(holderV.elementAt(checkIndex));
            checkIndex += 1;
            if (checkIndex == holderV.size()) {
                break;
            }
        }
        regionsV.add(holderV);
        if (regionsV.size() > 0x0000FFFF) {
            // unsigned short max, must match getShortEncodedRegionIndexImage()
            throw new ImageException("Error: image segmentation contains more than " + (0x0000FFFF) + " distinct regions.");
        }
        int regionSize = 0;
        byte pixelCheck = regionImagePixelV.elementAt(holderV.elementAt(0));
        for (int j = 0; j < holderV.size(); j++) {
            if (pixelCheck != regionImagePixelV.elementAt(holderV.elementAt(j))) {
                throw new ImageException("Linked regions have different image pixel values");
            }
            linkRegionMap[holderV.elementAt(j)] = regionsV.size() - 1;
            // System.out.print((j!=0?" ":"")+holderV.elementAt(j));
            regionSize += regionSizeArr[holderV.elementAt(j)];
            totalSize += regionSizeArr[holderV.elementAt(j)];
        }
        regionsSizeV.add(regionSize);
    // System.out.println();
    }
    if (totalSize != vcImage.getNumXYZ()) {
        throw new ImageException("Accumulated regions size does not equal image size");
    }
    // System.out.println("----------gather link distinct regions time "+((System.currentTimeMillis()-startTime)/1000.0));
    // startTime = System.currentTimeMillis();
    // //Create bitmasks of distinct regions
    // BitSet[] regionBitMaskBS = new BitSet[regionsV.size()];
    // for (int i = 0; i < regionBitMaskBS.length; i++) {
    // regionBitMaskBS[i] = new BitSet(regionPixels.length);
    // }
    // for (int i = 0; i < regionPixels.length; i++) {
    // regionBitMaskBS[linkRegionMap[regionPixels[i]]].set(i);
    // }
    // System.out.println("----------bitmask time "+((System.currentTimeMillis()-startTime)/1000.0));
    // startTime = System.currentTimeMillis();
    // Create RegionInfos
    regionInfos = new RegionInfo[regionsV.size()];
    for (int i = 0; i < regionsV.size(); i++) {
        regionInfos[i] = new RegionInfo(regionImagePixelV.elementAt(regionsV.elementAt(i).elementAt(0)) & 0x000000FF, regionsSizeV.elementAt(i), i, // regionBitMaskBS[i]
        null);
    }
    // System.out.println("----------regioninfo time "+((System.currentTimeMillis()-startTime)/1000.0));
    // startTime = System.currentTimeMillis();
    mapImageIndexToLinkRegion = regionPixels;
    mapLinkRegionToDistinctRegion = linkRegionMap;
    if (dimension != 0) {
        generateSurfaceCollection(regionsV.size(), // regionPixels,linkRegionMap,
        vcImage, xSurfElements, ySurfElements, zSurfElements, dimension, extent, origin);
    }
    if (surfaceCollection != null && debug_bCheckPolygonQuality) {
        verifyQuadVertexOrdering(debug_maxQuadAngle);
    }
    // Taubin smoothing
    if (surfaceCollection != null && filterCutoffFrequency < RegionImage.NO_SMOOTHING) {
        TaubinSmoothing taubinSmoothing = new TaubinSmoothingWrong();
        TaubinSmoothingSpecification taubinSpec = TaubinSmoothingSpecification.getInstance(filterCutoffFrequency);
        taubinSmoothing.smooth(surfaceCollection, taubinSpec);
    }
    // startTime = System.currentTimeMillis();
    if (surfaceCollection != null && debug_bCheckPolygonQuality) {
        verifyQuadVertexOrdering(debug_maxQuadAngle);
    }
// System.out.println("Total Num Regions = "+regionsV.size());
// System.out.println("Total Size = "+totalSize);
}