vtkSurfaceFromVolume.C 11.2 KB
Newer Older
hrchilds's avatar
hrchilds committed
1 2
/*****************************************************************************
*
3
* Copyright (c) 2000 - 2012, Lawrence Livermore National Security, LLC
hrchilds's avatar
hrchilds committed
4
* Produced at the Lawrence Livermore National Laboratory
5
* LLNL-CODE-442911
hrchilds's avatar
hrchilds committed
6 7
* All rights reserved.
*
8
* This file is  part of VisIt. For  details, see https://visit.llnl.gov/.  The
hrchilds's avatar
hrchilds committed
9 10 11 12 13 14 15 16 17 18
* full copyright notice is contained in the file COPYRIGHT located at the root
* of the VisIt distribution or at http://www.llnl.gov/visit/copyright.html.
*
* Redistribution  and  use  in  source  and  binary  forms,  with  or  without
* modification, are permitted provided that the following conditions are met:
*
*  - Redistributions of  source code must  retain the above  copyright notice,
*    this list of conditions and the disclaimer below.
*  - Redistributions in binary form must reproduce the above copyright notice,
*    this  list of  conditions  and  the  disclaimer (as noted below)  in  the
19 20 21
*    documentation and/or other materials provided with the distribution.
*  - Neither the name of  the LLNS/LLNL nor the names of  its contributors may
*    be used to endorse or promote products derived from this software without
hrchilds's avatar
hrchilds committed
22 23 24 25 26
*    specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT  HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR  PURPOSE
27 28 29
* ARE  DISCLAIMED. IN  NO EVENT  SHALL LAWRENCE  LIVERMORE NATIONAL  SECURITY,
* LLC, THE  U.S.  DEPARTMENT OF  ENERGY  OR  CONTRIBUTORS BE  LIABLE  FOR  ANY
* DIRECT,  INDIRECT,   INCIDENTAL,   SPECIAL,   EXEMPLARY,  OR   CONSEQUENTIAL
hrchilds's avatar
hrchilds committed
30 31 32 33 34 35 36 37 38
* DAMAGES (INCLUDING, BUT NOT  LIMITED TO, PROCUREMENT OF  SUBSTITUTE GOODS OR
* SERVICES; LOSS OF  USE, DATA, OR PROFITS; OR  BUSINESS INTERRUPTION) HOWEVER
* CAUSED  AND  ON  ANY  THEORY  OF  LIABILITY,  WHETHER  IN  CONTRACT,  STRICT
* LIABILITY, OR TORT  (INCLUDING NEGLIGENCE OR OTHERWISE)  ARISING IN ANY  WAY
* OUT OF THE  USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
* DAMAGE.
*
*****************************************************************************/

hrchilds's avatar
hrchilds committed
39 40 41 42 43 44
// ************************************************************************* //
//                           vtkSurfaceFromVolume.C                          //
// ************************************************************************* //

#include <vtkSurfaceFromVolume.h>

hrchilds's avatar
hrchilds committed
45
#include <vtkCellArray.h>
hrchilds's avatar
hrchilds committed
46
#include <vtkCellData.h>
hrchilds's avatar
hrchilds committed
47
#include <vtkIdTypeArray.h>
hrchilds's avatar
hrchilds committed
48
#include <vtkIntArray.h>
hrchilds's avatar
hrchilds committed
49 50 51
#include <vtkPointData.h>
#include <vtkPolyData.h>

52
#include <vtkAccessors.h>
hrchilds's avatar
hrchilds committed
53 54 55 56 57 58

vtkSurfaceFromVolume::TriangleList::TriangleList()
{
    listSize = 4096;
    trianglesPerList = 1024;
 
59 60 61
    list = new vtkIdType*[listSize];
    list[0] = new vtkIdType[4*trianglesPerList];
    for (vtkIdType i = 1 ; i < listSize ; i++)
hrchilds's avatar
hrchilds committed
62 63 64 65 66 67
        list[i] = NULL;
 
    currentList = 0;
    currentTriangle = 0;
}
 
hrchilds's avatar
hrchilds committed
68

hrchilds's avatar
hrchilds committed
69 70 71 72 73 74 75 76 77 78 79 80 81 82
vtkSurfaceFromVolume::TriangleList::~TriangleList()
{
    for (int i = 0 ; i < listSize ; i++)
    {
        if (list[i] != NULL)
            delete [] list[i];
        else
            break;
    }
    delete [] list;
}
 
 
int
83
vtkSurfaceFromVolume::TriangleList::GetList(int listId, const vtkIdType *&outlist)
hrchilds's avatar
hrchilds committed
84 85 86 87 88 89 90 91 92 93 94 95 96
    const
{
    if (listId < 0 || listId > currentList)
    {
        outlist = NULL;
        return 0;
    }
 
    outlist = list[listId];
    return (listId == currentList ? currentTriangle : trianglesPerList);
}
 
 
97
vtkIdType
hrchilds's avatar
hrchilds committed
98 99 100 101 102 103
vtkSurfaceFromVolume::TriangleList::GetNumberOfLists(void) const
{
    return currentList+1;
}


104
vtkIdType
hrchilds's avatar
hrchilds committed
105 106
vtkSurfaceFromVolume::TriangleList::GetTotalNumberOfTriangles(void) const
{
107 108
    vtkIdType numFullLists = currentList;  // actually currentList-1+1
    vtkIdType numExtra = currentTriangle;  // again, currentTriangle-1+1
hrchilds's avatar
hrchilds committed
109 110 111 112 113
 
    return numFullLists*trianglesPerList + numExtra;
}
 
 
hrchilds's avatar
hrchilds committed
114 115 116 117 118 119
// ****************************************************************************
//  Modifications:
//
//    Hank Childs, Thu Oct 21 15:32:07 PDT 2004
//    Fix bug when we get low on memory.
//
hrchilds's avatar
hrchilds committed
120 121 122
//    Sean Ahern, Mon Mar  5 15:44:05 EST 2007
//    Fixed test for resizing list.  Initialized new entries.
//
hrchilds's avatar
hrchilds committed
123 124
// ****************************************************************************

hrchilds's avatar
hrchilds committed
125
void
126 127
vtkSurfaceFromVolume::TriangleList::AddTriangle(vtkIdType cellId, 
    vtkIdType v1, vtkIdType v2, vtkIdType v3)
hrchilds's avatar
hrchilds committed
128 129 130
{
    if (currentTriangle >= trianglesPerList)
    {
hrchilds's avatar
hrchilds committed
131
        if ((currentList+1) >= listSize)
hrchilds's avatar
hrchilds committed
132
        {
133 134
            vtkIdType **tmpList = new vtkIdType*[2*listSize];
            for (vtkIdType i = 0 ; i < listSize ; i++)
hrchilds's avatar
hrchilds committed
135
                tmpList[i] = list[i];
136
            for (vtkIdType i = listSize ; i < listSize*2 ; i++)
hrchilds's avatar
hrchilds committed
137 138
                tmpList[i] = NULL;

hrchilds's avatar
hrchilds committed
139 140 141 142 143 144
            listSize *= 2;
            delete [] list;
            list = tmpList;
        }
 
        currentList++;
145
        list[currentList] = new vtkIdType[4*trianglesPerList];
hrchilds's avatar
hrchilds committed
146 147 148
        currentTriangle = 0;
    }
 
149
    vtkIdType idx = 4*currentTriangle;
hrchilds's avatar
hrchilds committed
150 151 152 153 154 155 156
    list[currentList][idx+0] = cellId;
    list[currentList][idx+1] = v1;
    list[currentList][idx+2] = v2;
    list[currentList][idx+3] = v3;
    currentTriangle++;
}

hrchilds's avatar
hrchilds committed
157 158
// ****************************************************************************
//  Modications:
159 160 161 162 163
//    Brad Whitlock, Mon Mar 12 16:15:37 PDT 2012
//    I separated this code out from ConstructPolyData into its own function so 
//    I could parameterize the point access code using a "point getter"
//    object. This code gets inlined in the ConstructPolyData methods so we
//    can support multiple coordinate data types.
hrchilds's avatar
hrchilds committed
164
//
hrchilds's avatar
hrchilds committed
165
// ****************************************************************************
hrchilds's avatar
hrchilds committed
166

167 168 169 170 171 172
template <typename PointGetter>
inline void
ConstructPolyDataHelper(vtkPointData *inPD, vtkCellData *inCD,
     vtkPolyData *output, const vtkSurfaceFromVolume::PointList &pt_list,
     const vtkSurfaceFromVolume::TriangleList &tris,
     int dataType, const PointGetter &pointGetter)
hrchilds's avatar
hrchilds committed
173
{
174
    vtkIdType   i, j;
hrchilds's avatar
hrchilds committed
175 176 177 178

    vtkPointData *outPD = output->GetPointData();
    vtkCellData  *outCD = output->GetCellData();

hrchilds's avatar
hrchilds committed
179 180 181 182
    vtkIntArray *newOrigNodes = NULL;
    vtkIntArray *origNodes = vtkIntArray::SafeDownCast(
              inPD->GetArray("avtOriginalNodeNumbers"));

hrchilds's avatar
hrchilds committed
183 184 185
    //
    // Set up the output points and its point data.
    //
186 187
    vtkPoints *outPts = vtkPoints::New(dataType);
    vtkIdType nOutPts = pt_list.GetTotalNumberOfPoints();
hrchilds's avatar
hrchilds committed
188 189
    outPts->SetNumberOfPoints(nOutPts);
    outPD->CopyAllocate(inPD, nOutPts);
hrchilds's avatar
hrchilds committed
190 191 192 193 194 195 196
    if (origNodes != NULL)
    {
        newOrigNodes = vtkIntArray::New();
        newOrigNodes->SetNumberOfComponents(origNodes->GetNumberOfComponents());
        newOrigNodes->SetNumberOfTuples(nOutPts);
        newOrigNodes->SetName(origNodes->GetName());
    }
197 198
    vtkIdType nLists = pt_list.GetNumberOfLists();
    vtkIdType ptIdx = 0;
hrchilds's avatar
hrchilds committed
199 200
    for (i = 0 ; i < nLists ; i++)
    {
201 202
        const vtkDataSetFromVolume::PointEntry *pe_list = NULL;
        vtkIdType nPts = pt_list.GetList(i, pe_list);
hrchilds's avatar
hrchilds committed
203 204
        for (j = 0 ; j < nPts ; j++)
        {
205 206 207 208 209 210 211 212 213
            const vtkDataSetFromVolume::PointEntry &pe = pe_list[j];
            double pt[3], pt1[3], pt2[3];
            pointGetter.GetPoint(pe.ptIds[0], pt1);
            pointGetter.GetPoint(pe.ptIds[1], pt2);
            double p  = pe.percent;
            double bp = 1. - p;
            pt[0] = pt1[0]*p + pt2[0]*bp;
            pt[1] = pt1[1]*p + pt2[1]*bp;
            pt[2] = pt1[2]*p + pt2[2]*bp;
hrchilds's avatar
hrchilds committed
214
            outPts->SetPoint(ptIdx, pt);
hrchilds's avatar
hrchilds committed
215
            outPD->InterpolateEdge(inPD, ptIdx, pe.ptIds[0], pe.ptIds[1], bp);
hrchilds's avatar
hrchilds committed
216 217
            if (newOrigNodes)
            {
218
                vtkIdType id = (bp <= 0.5 ? pe.ptIds[0] : pe.ptIds[1]);
hrchilds's avatar
hrchilds committed
219 220
                newOrigNodes->SetTuple(ptIdx, origNodes->GetTuple(id));
            }
hrchilds's avatar
hrchilds committed
221 222 223 224 225
            ptIdx++;
        }
    }
    output->SetPoints(outPts);
    outPts->Delete();
hrchilds's avatar
hrchilds committed
226 227 228 229 230 231 232
    if (newOrigNodes)
    {
        // AddArray will overwrite an already existing array with 
        // the same name, exactly what we want here.
        outPD->AddArray(newOrigNodes);
        newOrigNodes->Delete();
    }
hrchilds's avatar
hrchilds committed
233 234 235 236

    //
    // Now set up the triangles and the cell data.
    //
237
    vtkIdType ntris = tris.GetTotalNumberOfTriangles();
hrchilds's avatar
hrchilds committed
238 239 240 241
    vtkIdTypeArray *nlist = vtkIdTypeArray::New();
    nlist->SetNumberOfValues(3*ntris + ntris);
    vtkIdType *nl = nlist->GetPointer(0);

hrchilds's avatar
hrchilds committed
242
    outCD->CopyAllocate(inCD, ntris);
243 244
    vtkIdType cellId = 0;
    vtkIdType nlists = tris.GetNumberOfLists();
hrchilds's avatar
hrchilds committed
245 246
    for (i = 0 ; i < nlists ; i++)
    {
247 248
        const vtkIdType *list;
        vtkIdType listSize = tris.GetList(i, list);
hrchilds's avatar
hrchilds committed
249 250 251
        for (j = 0 ; j < listSize ; j++)
        {
            outCD->CopyData(inCD, list[0], cellId);
hrchilds's avatar
hrchilds committed
252 253 254 255
            *nl++ = 3;
            *nl++ = list[1];
            *nl++ = list[2];
            *nl++ = list[3];
hrchilds's avatar
hrchilds committed
256 257 258 259
            list += 4;
            cellId++;
        }
    }
hrchilds's avatar
hrchilds committed
260 261 262 263 264 265
    vtkCellArray *cells = vtkCellArray::New();
    cells->SetCells(ntris, nlist);
    nlist->Delete();

    output->SetPolys(cells);
    cells->Delete();
hrchilds's avatar
hrchilds committed
266 267
}

268 269 270 271 272 273 274 275 276 277 278
// ****************************************************************************
//  Modications:
//
//    Hank Childs, Fri Jan 30 08:50:44 PST 2004
//    Speed up construction by doing pointer arithmetic.
//
//    Kathleen Bonnell, Fri Nov 12 16:43:11 PST 2004 
//    Don't interoplate avtOriginalNodeNumbers, instead use the closest
//    point to the slice point.
//
// ****************************************************************************
hrchilds's avatar
hrchilds committed
279

280 281 282
void
vtkSurfaceFromVolume::ConstructPolyData(vtkPointData *inPD, vtkCellData *inCD,
                                        vtkPolyData *output, vtkPoints *pts)
hrchilds's avatar
hrchilds committed
283
{
284 285 286 287 288 289
    if(pts->GetDataType() == VTK_FLOAT)
        ConstructPolyDataHelper(inPD, inCD, output, this->pt_list, this->tris, pts->GetDataType(), vtkPointAccessor<float>(pts));
    else if(pts->GetDataType() == VTK_FLOAT)
        ConstructPolyDataHelper(inPD, inCD, output, this->pt_list, this->tris, pts->GetDataType(), vtkPointAccessor<double>(pts));
    else
        ConstructPolyDataHelper(inPD, inCD, output, this->pt_list, this->tris, pts->GetDataType(), vtkGeneralPointAccessor(pts));
hrchilds's avatar
hrchilds committed
290 291
}

hrchilds's avatar
hrchilds committed
292 293 294 295 296 297
// ****************************************************************************
//  Modications:
//
//    Hank Childs, Fri Jan 30 08:50:44 PST 2004
//    Speed up construction by doing pointer arithmetic.
//
hrchilds's avatar
hrchilds committed
298 299 300 301
//    Kathleen Bonnell, Tue Nov 23 14:07:10 PST 2004 
//    Don't interoplate avtOriginalNodeNumbers, instead use the closest
//    point to the slice point.
//
302 303 304 305
//    Brad Whitlock, Mon Mar 12 16:09:19 PDT 2012
//    I moved indexing code into "point accessor" classes and added support
//    for different coordinate types.
//
hrchilds's avatar
hrchilds committed
306 307
// ****************************************************************************

hrchilds's avatar
hrchilds committed
308 309
void
vtkSurfaceFromVolume::ConstructPolyData(vtkPointData *inPD, vtkCellData *inCD,
310 311
    vtkPolyData *output, int *dims, 
    vtkDataArray *X, vtkDataArray *Y, vtkDataArray *Z)
hrchilds's avatar
hrchilds committed
312
{
313 314 315 316 317 318 319 320 321 322
    int tx = X->GetDataType();
    int ty = Y->GetDataType();
    int tz = Z->GetDataType();
    bool same = tx == ty && ty == tz;
    if(same && tx == VTK_FLOAT)
        ConstructPolyDataHelper(inPD, inCD, output, this->pt_list, this->tris, tx, vtkRectPointAccessor<float>(dims, X, Y, Z));
    else if(same && tx == VTK_DOUBLE)
        ConstructPolyDataHelper(inPD, inCD, output, this->pt_list, this->tris, tx, vtkRectPointAccessor<double>(dims, X, Y, Z));
    else
        ConstructPolyDataHelper(inPD, inCD, output, this->pt_list, this->tris, tx, vtkGeneralRectPointAccessor(dims, X, Y, Z));
hrchilds's avatar
hrchilds committed
323
}