1 /*
2  * Medical Image Registration ToolKit (MIRTK)
3  *
4  * Copyright 2013-2015 Imperial College London
5  * Copyright 2013-2015 Andreas Schuh
6  *
7  * Licensed under the Apache License, Version 2.0 (the "License");
8  * you may not use this file except in compliance with the License.
9  * You may obtain a copy of the License at
10  *
11  *     http://www.apache.org/licenses/LICENSE-2.0
12  *
13  * Unless required by applicable law or agreed to in writing, software
14  * distributed under the License is distributed on an "AS IS" BASIS,
15  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16  * See the License for the specific language governing permissions and
17  * limitations under the License.
18  */
19 
20 #include "mirtk/Common.h"
21 #include "mirtk/Options.h"
22 
23 #include "mirtk/PointSetIO.h"
24 #include "mirtk/PointSetUtils.h"
25 
26 #include "vtkNew.h"
27 #include "vtkSmartPointer.h"
28 #include "vtkPointSet.h"
29 #include "vtkPointData.h"
30 #include "vtkCellData.h"
31 #include "vtkDataArray.h"
32 #include "vtkFloatArray.h"
33 #include "vtkDataSetAttributes.h"
34 #include "vtkPointDataToCellData.h"
35 #include "vtkCellDataToPointData.h"
36 
37 using namespace mirtk;
38 
39 
40 // =============================================================================
41 // Help
42 // =============================================================================
43 
44 // -----------------------------------------------------------------------------
PrintHelp(const char * name)45 void PrintHelp(const char *name)
46 {
47   cout << endl;
48   cout << "Usage: " << name << " <source> <target> [<output>] [options]" << endl;
49   cout << endl;
50   cout << "Description:" << endl;
51   cout << "  Copies point and/or cell data from a source point set to a target point set" << endl;
52   cout << "  and writes the resulting amended point set to the specified output file." << endl;
53   cout << "  When no separate output file is specified, the target point set is overwritten." << endl;
54   cout << "  This command can also convert from point data to cell data and vice versa." << endl;
55   cout << endl;
56   cout << "  If the point sets have differing number of points or cells, respectively," << endl;
57   cout << "  zero entries are either added to the target arrays or only the first n tuples" << endl;
58   cout << "  of the source arrays are copied. In case of the -points option, the remaining" << endl;
59   cout << "  target points are kept unchanged." << endl;
60   cout << endl;
61   cout << "  If the <attr> argument is given to any of the option below, the copied data" << endl;
62   cout << "  array is set as the new \"scalars\", \"vectors\", \"tcoords\", or \"normals\"" << endl;
63   cout << "  attribute, respectively, of the output point set." << endl;
64   cout << endl;
65   cout << "  When no options are given, the :option:`-scalars` are copied over." << endl;
66   cout << endl;
67   cout << "Arguments:" << endl;
68   cout << "  source   Point set from which to copy the point/cell data." << endl;
69   #if MIRTK_IO_WITH_GIFTI
70   cout << "           Can also be a GIFTI file with only point data arrays." << endl;
71   cout << "           In this case, the coordinates and topology of the target" << endl;
72   cout << "           surface are used to define the source surface." << endl;
73   #endif
74   cout << "  target   Point set to which to add the point/cell data." << endl;
75   cout << "  output   Output point set with copied point/cell data." << endl;
76   cout << "           If not specified, the target file is overwritten." << endl;
77   cout << endl;
78   cout << "Optional arguments:" << endl;
79   cout << "  -scalars" << endl;
80   cout << "      Copy \"scalars\" data set attribute. (default)" << endl;
81   cout << endl;
82   cout << "  -vectors" << endl;
83   cout << "      Copy \"vectors\" data set attribute." << endl;
84   cout << endl;
85   cout << "  -vectors-as-points" << endl;
86   cout << "      Set points of target to \"vectors\" data set attribute of source." << endl;
87   cout << endl;
88   cout << "  -normals" << endl;
89   cout << "      Copy \"normals\" data set attribute." << endl;
90   cout << endl;
91   cout << "  -tcoords" << endl;
92   cout << "      Copy \"tcoords\" data set attribute." << endl;
93   cout << endl;
94   cout << "  -tcoords-as-points" << endl;
95   cout << "      Set points of target to \"tcoords\" data set attribute of source." << endl;
96   cout << endl;
97   cout << "  -points [<target_name> [<attr>]]" << endl;
98   cout << "      Add points of source as point data of target." << endl;
99   cout << "      When no <target_name> (and <attr>) specified, the points of the" << endl;
100   cout << "      target point set are replaced by those of the source point set." << endl;
101   cout << "      If the target point set has more points than the source point set," << endl;
102   cout << "      the remaining points of the target point set are unchanged." << endl;
103   cout << endl;
104   cout << "  -pointdata <name|index> [<target_name> [<attr>]]" << endl;
105   cout << "      Name or index of source point data to copy." << endl;
106   cout << endl;
107   cout << "  -pointdata-as-points <name|index>" << endl;
108   cout << "      Replaces the points of the target point set by the first three" << endl;
109   cout << "      components of the specified source point data array." << endl;
110   cout << endl;
111   cout << "  -pointdata-as-celldata <name|index>" << endl;
112   cout << "      Converts the specified point data array of the source point set to cell" << endl;
113   cout << "      data and adds a corresponding cell data array to the target point set." << endl;
114   cout << endl;
115   cout << "  -pointmask <name>" << endl;
116   cout << "      Add point data array which indicates copied point data." << endl;
117   cout << endl;
118   cout << "  -celldata  <name|index> [<target_name> [<attr>]]" << endl;
119   cout << "      Name or index of source cell data to copy." << endl;
120   cout << endl;
121   cout << "  -celldata-as-pointdata <name|index>" << endl;
122   cout << "      Converts the specified cell data array of the source point set to point" << endl;
123   cout << "      data and adds a corresponding point data array to the target point set." << endl;
124   cout << endl;
125   cout << "  -cellmask <name>" << endl;
126   cout << "      Add cell data array which indicates copied cell data." << endl;
127   cout << endl;
128   cout << "  -case-[in]sensitive" << endl;
129   cout << "      Lookup source arrays by case [in]sensitive name. (default: insensitive)" << endl;
130   cout << endl;
131   cout << "  -majority" << endl;
132   cout << "      Use majority vote for resampling categorical data. (default)" << endl;
133   cout << endl;
134   cout << "  -unanimous" << endl;
135   cout << "      Use unanimous vote for resampling categorical data." << endl;
136   PrintCommonOptions(cout);
137   cout << endl;
138 }
139 
140 // =============================================================================
141 // Auxiliaries
142 // =============================================================================
143 
144 // -----------------------------------------------------------------------------
145 using AttributeType                = vtkDataSetAttributes::AttributeTypes;
146 const AttributeType NUM_ATTRIBUTES = vtkDataSetAttributes::NUM_ATTRIBUTES;
147 
148 // -----------------------------------------------------------------------------
149 struct ArrayInfo
150 {
151   int           _SourceIndex;
152   const char   *_SourceName;
153   AttributeType _SourceAttribute;
154   int           _TargetIndex;
155   const char   *_TargetName;
156   AttributeType _TargetAttribute;
157   int           _TargetDataType;
158   bool          _PointCellConversion;
159 
ArrayInfoArrayInfo160   ArrayInfo()
161   :
162     _SourceIndex(-1),
163     _SourceName(nullptr),
164     _SourceAttribute(NUM_ATTRIBUTES),
165     _TargetIndex(-1),
166     _TargetName(nullptr),
167     _TargetAttribute(NUM_ATTRIBUTES),
168     _TargetDataType(VTK_VOID),
169     _PointCellConversion(false)
170   {}
171 };
172 
173 // -----------------------------------------------------------------------------
174 /// Get source point data array
GetSourceArray(const char * type,vtkDataSetAttributes * attr,const ArrayInfo & info,bool case_sensitive)175 vtkDataArray *GetSourceArray(const char *type, vtkDataSetAttributes *attr, const ArrayInfo &info, bool case_sensitive)
176 {
177   vtkDataArray *array;
178   if (info._SourceAttribute < NUM_ATTRIBUTES) {
179     array = attr->GetAttribute(info._SourceAttribute);
180   } else if (info._SourceName) {
181     if (case_sensitive) {
182       array = attr->GetArray(info._SourceName);
183     } else {
184       array = GetArrayByCaseInsensitiveName(attr, info._SourceName);
185     }
186     if (array == nullptr) {
187       FatalError("Source has no " << type << " data array named " << info._SourceName);
188     }
189   } else {
190     if (info._SourceIndex < 0 || info._SourceIndex >= attr->GetNumberOfArrays()) {
191       FatalError("Source has no " << type << " data array with index " << info._SourceIndex);
192     }
193     array = attr->GetArray(info._SourceIndex);
194   }
195   return array;
196 }
197 
198 // -----------------------------------------------------------------------------
199 /// Parse VTK data type string
ParseVtkDataType(const char * arg)200 int ParseVtkDataType(const char *arg)
201 {
202   const string lstr = ToLower(Trim(arg));
203   if (lstr == "char")   return VTK_CHAR;
204   if (lstr == "uchar")  return VTK_UNSIGNED_CHAR;
205   if (lstr == "binary") return VTK_UNSIGNED_CHAR;
206   if (lstr == "short")  return VTK_SHORT;
207   if (lstr == "grey")   return VTK_SHORT;
208   if (lstr == "ushort") return VTK_UNSIGNED_SHORT;
209   if (lstr == "int")    return VTK_INT;
210   if (lstr == "uint")   return VTK_UNSIGNED_INT;
211   if (lstr == "long")   return VTK_LONG_LONG;
212   if (lstr == "ulong")  return VTK_UNSIGNED_LONG_LONG;
213   if (lstr == "float")  return VTK_FLOAT;
214   if (lstr == "double") return VTK_DOUBLE;
215   if (lstr == "real") {
216     #if MIRTK_USE_FLOAT_BY_DEFAULT
217       return VTK_FLOAT;
218     #else
219       return VTK_DOUBLE;
220     #endif
221   }
222   FatalError("Unknown data type: " << arg);
223 }
224 
225 /// Enumeration of categorical data interpolation methods
226 enum LabelInterpolationMode
227 {
228   MajorityVote,
229   UnanimousVote
230 };
231 
232 // -----------------------------------------------------------------------------
233 /// Convert point data labels to cell data
234 vtkSmartPointer<vtkDataArray>
ConvertPointToCellLabels(vtkPointSet * pset,vtkDataArray * labels,LabelInterpolationMode mode=MajorityVote)235 ConvertPointToCellLabels(vtkPointSet *pset, vtkDataArray *labels,
236                          LabelInterpolationMode mode = MajorityVote)
237 {
238   vtkSmartPointer<vtkDataArray> output;
239   output.TakeReference(labels->NewInstance());
240   output->SetName(labels->GetName());
241   for (int j = 0; j < labels->GetNumberOfComponents(); ++j) {
242     output->SetComponentName(j, labels->GetComponentName(j));
243   }
244   output->SetNumberOfComponents(labels->GetNumberOfComponents());
245   output->SetNumberOfTuples(pset->GetNumberOfCells());
246 
247   switch (mode) {
248     case MajorityVote: {
249       int    count;
250       double label;
251       OrderedMap<double, int> hist;
252       OrderedMap<double, int>::iterator bin;
253       vtkNew<vtkIdList> ptIds;
254       for (vtkIdType cellId = 0; cellId < pset->GetNumberOfCells(); ++cellId) {
255         pset->GetCellPoints(cellId, ptIds.GetPointer());
256         for (int j = 0; j < labels->GetNumberOfComponents(); ++j) {
257           hist.clear();
258           for (vtkIdType i = 0; i < ptIds->GetNumberOfIds(); ++i) {
259             label = labels->GetComponent(ptIds->GetId(i), j);
260             bin = hist.find(label);
261             if (bin == hist.end()) hist[label] = 1;
262             else bin->second += 1;
263           }
264           label = 0., count = 0;
265           for (bin = hist.begin(); bin != hist.end(); ++bin) {
266             if (bin->second > count) {
267               label = bin->first;
268               count = bin->second;
269             }
270           }
271           output->SetComponent(cellId, j, label);
272         }
273       }
274     } break;
275     case UnanimousVote: {
276       const double invalid = NaN;
277       double       label;
278       vtkNew<vtkIdList> ptIds;
279       for (vtkIdType cellId = 0; cellId < pset->GetNumberOfCells(); ++cellId) {
280         pset->GetCellPoints(cellId, ptIds.GetPointer());
281         for (int j = 0; j < labels->GetNumberOfComponents(); ++j) {
282           if (ptIds->GetNumberOfIds() == 0) {
283             label = invalid;
284           } else {
285             label = labels->GetComponent(ptIds->GetId(0), j);
286             for (vtkIdType i = 1; i < ptIds->GetNumberOfIds(); ++i) {
287               if (labels->GetComponent(ptIds->GetId(i), j) != label) {
288                 label = invalid;
289                 break;
290               }
291             }
292           }
293           output->SetComponent(cellId, j, label);
294         }
295       }
296     } break;
297   }
298 
299   return output;
300 }
301 
302 // -----------------------------------------------------------------------------
303 /// Convert cell data labels to point data
304 vtkSmartPointer<vtkDataArray>
ConvertCellToPointLabels(vtkPointSet * pset,vtkDataArray * labels,LabelInterpolationMode mode=MajorityVote)305 ConvertCellToPointLabels(vtkPointSet *pset, vtkDataArray *labels,
306                          LabelInterpolationMode mode = MajorityVote)
307 {
308   vtkSmartPointer<vtkDataArray> output;
309   output.TakeReference(labels->NewInstance());
310   output->SetName(labels->GetName());
311   for (int j = 0; j < labels->GetNumberOfComponents(); ++j) {
312     output->SetComponentName(j, labels->GetComponentName(j));
313   }
314   output->SetNumberOfComponents(labels->GetNumberOfComponents());
315   output->SetNumberOfTuples(pset->GetNumberOfPoints());
316 
317   switch (mode) {
318     case MajorityVote: {
319       int    count;
320       double label;
321       OrderedMap<double, int> hist;
322       OrderedMap<double, int>::iterator bin;
323       vtkNew<vtkIdList> cellIds;
324       for (vtkIdType ptId = 0; ptId < pset->GetNumberOfPoints(); ++ptId) {
325         pset->GetPointCells(ptId, cellIds.GetPointer());
326         for (int j = 0; j < labels->GetNumberOfComponents(); ++j) {
327           hist.clear();
328           for (vtkIdType i = 0; i < cellIds->GetNumberOfIds(); ++i) {
329             label = labels->GetComponent(cellIds->GetId(i), j);
330             bin = hist.find(label);
331             if (bin == hist.end()) hist[label] = 1;
332             else bin->second += 1;
333           }
334           label = 0., count = 0;
335           for (bin = hist.begin(); bin != hist.end(); ++bin) {
336             if (bin->second > count) {
337               label = bin->first;
338               count = bin->second;
339             }
340           }
341           output->SetComponent(ptId, j, label);
342         }
343       }
344     } break;
345     case UnanimousVote: {
346       double       label;
347       const double invalid = NaN;
348       vtkNew<vtkIdList> cellIds;
349       for (vtkIdType ptId = 0; ptId < pset->GetNumberOfPoints(); ++ptId) {
350         pset->GetPointCells(ptId, cellIds.GetPointer());
351         for (int j = 0; j < labels->GetNumberOfComponents(); ++j) {
352           if (cellIds->GetNumberOfIds() == 0) {
353             label = invalid;
354           } else {
355             label = labels->GetComponent(cellIds->GetId(0), j);
356             for (vtkIdType i = 1; i < cellIds->GetNumberOfIds(); ++i) {
357               if (labels->GetComponent(cellIds->GetId(i), j) != label) {
358                 label = invalid;
359                 break;
360               }
361             }
362           }
363           output->SetComponent(ptId, j, label);
364         }
365       }
366     } break;
367   }
368 
369   return output;
370 }
371 
372 // -----------------------------------------------------------------------------
373 /// Copy tuples from the input data array
Copy(vtkDataArray * array,vtkIdType n,int dtype=VTK_VOID)374 vtkSmartPointer<vtkDataArray> Copy(vtkDataArray *array, vtkIdType n, int dtype = VTK_VOID)
375 {
376   vtkSmartPointer<vtkDataArray> copy;
377   if (dtype == VTK_VOID) {
378     copy.TakeReference(array->NewInstance());
379   } else {
380     copy = NewVtkDataArray(dtype);
381   }
382   copy->SetName(array->GetName());
383   copy->SetNumberOfComponents(array->GetNumberOfComponents());
384   copy->SetNumberOfTuples(n);
385   for (int j = 0; j < array->GetNumberOfComponents(); ++j) {
386     copy->SetComponentName(j, array->GetComponentName(j));
387   }
388   vtkIdType ptId = 0;
389   while (ptId < n && ptId < array->GetNumberOfTuples()) {
390     copy->SetTuple(ptId, array->GetTuple(ptId));
391     ++ptId;
392   }
393   while (ptId < n) {
394     for (int j = 0; j < copy->GetNumberOfComponents(); ++j) {
395       copy->SetComponent(ptId, j, .0);
396     }
397     ++ptId;
398   }
399   return copy;
400 }
401 
402 // -----------------------------------------------------------------------------
403 /// Add data array to data set attributes
AddArray(vtkDataSetAttributes * data,vtkDataArray * array,AttributeType type)404 void AddArray(vtkDataSetAttributes *data, vtkDataArray *array, AttributeType type)
405 {
406   switch (type) {
407     case AttributeType::SCALARS:     data->SetScalars(array); break;
408     case AttributeType::VECTORS:     data->SetVectors(array); break;
409     case AttributeType::NORMALS:     data->SetNormals(array); break;
410     case AttributeType::TCOORDS:     data->SetTCoords(array); break;
411     case AttributeType::TENSORS:     data->SetTensors(array); break;
412     case AttributeType::GLOBALIDS:   data->SetGlobalIds(array); break;
413     case AttributeType::PEDIGREEIDS: data->SetPedigreeIds(array); break;
414     default: data->AddArray(array); break;
415   }
416 }
417 
418 // =============================================================================
419 // Main
420 // =============================================================================
421 
422 // -----------------------------------------------------------------------------
main(int argc,char ** argv)423 int main(int argc, char **argv)
424 {
425   vtkSmartPointer<vtkDataArray> array;
426   vtkSmartPointer<vtkDataArray> copy;
427   vtkDataSetAttributes         *sourceAttr, *targetAttr;
428   vtkIdType                     ntuples;
429 
430   // Positional arguments
431   REQUIRES_POSARGS(2);
432 
433   if (NUM_POSARGS > 3) {
434     PrintHelp(EXECNAME);
435     exit(1);
436   }
437 
438   const char *source_name = POSARG(1);
439   const char *target_name = POSARG(2);
440   const char *output_name = (NUM_POSARGS == 3 ? POSARG(3) : target_name);
441 
442   // Optional arguments
443   Array<ArrayInfo> pd, cd;
444   const char *point_mask_name = nullptr;
445   const char *cell_mask_name  = nullptr;
446   bool        case_sensitive  = false;
447   LabelInterpolationMode label_mode = MajorityVote;
448 
449   for (ALL_OPTIONS) {
450     if (OPTION("-points")) {
451       ArrayInfo info;
452       info._SourceIndex = -2;
453       if (HAS_ARGUMENT) info._TargetName = ARGUMENT;
454       else              info._TargetName = NULL;
455       if (HAS_ARGUMENT) PARSE_ARGUMENT(info._TargetAttribute);
456       pd.push_back(info);
457     }
458     else if (OPTION("-scalars") || OPTION("-vectors") || OPTION("-normals") || OPTION("-tcoords")) {
459       ArrayInfo info;
460       FromString(OPTNAME + 1, info._SourceAttribute);
461       pd.push_back(info);
462       cd.push_back(info);
463     }
464     else if (OPTION("-vectors-as-points") || OPTION("-tcoords-as-points")) {
465       ArrayInfo info;
466       FromString(OPTNAME + 1, info._SourceAttribute);
467       info._TargetIndex = -2;
468       pd.push_back(info);
469     }
470     else if (OPTION("-pd") || OPTION("-pointdata") || OPTION("-pd-as-cd") || OPTION("-pointdata-as-celldata") ||
471              OPTION("-cd") || OPTION("-celldata")  || OPTION("-cd-as-pd") || OPTION("-celldata-as-pointdata")) {
472       ArrayInfo info;
473       if (strcmp(OPTNAME, "-pointdata-as-celldata") == 0 || strcmp(OPTNAME, "-pd-as-cd") == 0 ||
474           strcmp(OPTNAME, "-celldata-as-pointdata") == 0 || strcmp(OPTNAME, "-cd-as-pd") == 0) {
475         info._PointCellConversion = true;
476       } else {
477         info._PointCellConversion = false;
478       }
479       info._SourceName = ARGUMENT;
480       if (FromString(info._SourceName, info._SourceIndex)) {
481         info._SourceName = NULL;
482       } else {
483         info._SourceIndex = -1;
484       }
485       if (HAS_ARGUMENT) {
486         info._TargetName = ARGUMENT;
487       }
488       if (HAS_ARGUMENT) PARSE_ARGUMENT(info._TargetAttribute);
489       if (HAS_ARGUMENT) info._TargetDataType = ParseVtkDataType(ARGUMENT);
490       if (strncmp(OPTNAME, "-pointdata", 10) == 0 || strncmp(OPTNAME, "-pd", 3) == 0) {
491         pd.push_back(info);
492       } else {
493         cd.push_back(info);
494       }
495     }
496     else if (OPTION("-pointdata-as-points") || OPTION("-pd-as-points")) {
497       ArrayInfo info;
498       info._SourceName = ARGUMENT;
499       if (FromString(info._SourceName, info._SourceIndex)) {
500         info._SourceName = NULL;
501       } else {
502         info._SourceIndex = -1;
503       }
504       info._TargetIndex = -2;
505       pd.push_back(info);
506     }
507     else if (OPTION("-majority"))  label_mode = MajorityVote;
508     else if (OPTION("-unanimous")) label_mode = UnanimousVote;
509     else if (OPTION("-pointmask")) point_mask_name = ARGUMENT;
510     else if (OPTION("-cellmask"))  cell_mask_name  = ARGUMENT;
511     else if (OPTION("-case-sensitive"))   case_sensitive = true;
512     else if (OPTION("-case-insensitive")) case_sensitive = false;
513     else HANDLE_COMMON_OR_UNKNOWN_OPTION();
514   }
515 
516   // Read input data sets
517   vtkSmartPointer<vtkPointSet> source, target;
518   target = ReadPointSet(target_name);
519   #if MIRTK_IO_WITH_GIFTI
520     const string source_ext = Extension(source_name, EXT_Last);
521     if (source_ext == ".gii") {
522       source = ReadGIFTI(source_name, vtkPolyData::SafeDownCast(target));
523     } else if (source_ext == ".csv" || source_ext == ".tsv" || source_ext == ".txt") {
524       char sep = ',';
525       if (source_ext == ".tsv") sep = '\t';
526       source = ReadPointSetTable(source_name, sep, target);
527     }
528   #endif
529   if (source == nullptr) {
530     source = ReadPointSet(source_name);
531   }
532 
533   const vtkIdType npoints = target->GetNumberOfPoints();
534   const vtkIdType ncells  = target->GetNumberOfCells();
535 
536   vtkPointData * const sourcePD = source->GetPointData();
537   vtkCellData  * const sourceCD = source->GetCellData();
538 
539   vtkPointData * const targetPD = target->GetPointData();
540   vtkCellData  * const targetCD = target->GetCellData();
541 
542   // Copy all point set attributes by default
543   if (pd.empty() && cd.empty()) {
544     int attr;
545     for (int i = 0; i < sourcePD->GetNumberOfArrays(); ++i) {
546       attr = sourcePD->IsArrayAnAttribute(i);
547       if (attr < 0) attr = NUM_ATTRIBUTES;
548       ArrayInfo info;
549       info._SourceIndex     = i;
550       info._TargetAttribute = static_cast<AttributeType>(attr);
551       pd.push_back(info);
552     }
553     for (int i = 0; i < sourceCD->GetNumberOfArrays(); ++i) {
554       attr = sourceCD->IsArrayAnAttribute(i);
555       if (attr < 0) attr = NUM_ATTRIBUTES;
556       ArrayInfo info;
557       info._SourceIndex     = i;
558       info._TargetAttribute = static_cast<AttributeType>(attr);
559       cd.push_back(info);
560     }
561   }
562 
563   // Convert point data to cell data if necessary
564   vtkSmartPointer<vtkCellData> pd_as_cd;
565   for (size_t i = 0; i < pd.size(); ++i) {
566     if (pd[i]._PointCellConversion) {
567       vtkNew<vtkPointDataToCellData> pd_to_cd;
568       SetVTKInput(pd_to_cd, source);
569       pd_to_cd->PassPointDataOff();
570       pd_to_cd->Update();
571       pd_as_cd = pd_to_cd->GetOutput()->GetCellData();
572       break;
573     }
574   }
575 
576   // Convert cell data to point data if necessary
577   vtkSmartPointer<vtkPointData> cd_as_pd;
578   for (size_t i = 0; i < cd.size(); ++i) {
579     if (cd[i]._PointCellConversion) {
580       vtkNew<vtkCellDataToPointData> cd_to_pd;
581       SetVTKInput(cd_to_pd, source);
582       cd_to_pd->PassCellDataOff();
583       cd_to_pd->Update();
584       cd_as_pd = cd_to_pd->GetOutput()->GetPointData();
585       break;
586     }
587   }
588 
589   // Add point data
590   for (size_t i = 0; i < pd.size(); ++i) {
591     if (pd[i]._SourceIndex == -2) {
592       vtkIdType end = min(npoints, source->GetNumberOfPoints());
593       if (pd[i]._TargetName) {
594         copy = vtkSmartPointer<vtkFloatArray>::New();
595         copy->SetNumberOfComponents(3);
596         copy->SetNumberOfTuples(npoints);
597         for (vtkIdType ptId = 0; ptId < end; ++ptId) {
598           copy->SetTuple(ptId, source->GetPoint(ptId));
599         }
600         const double zero[3] = {.0, .0, .0};
601         for (vtkIdType ptId = end; ptId < npoints; ++ptId) {
602           copy->SetTuple(ptId, zero);
603         }
604       } else {
605         vtkPoints *points = target->GetPoints();
606         for (vtkIdType ptId = 0; ptId < end; ++ptId) {
607           points->SetPoint(ptId, source->GetPoint(ptId));
608         }
609         copy = nullptr;
610       }
611     } else {
612       if (pd[i]._PointCellConversion && pd[i]._TargetIndex != -2) {
613         sourceAttr = pd_as_cd;
614         targetAttr = targetCD;
615         ntuples    = ncells;
616       } else {
617         sourceAttr = sourcePD;
618         targetAttr = targetPD;
619         ntuples    = npoints;
620       }
621       array = GetSourceArray("point", sourceAttr, pd[i], case_sensitive);
622       if (pd[i]._TargetIndex == -2) {
623         double p[3] = {.0};
624         vtkPoints *points = target->GetPoints();
625         vtkIdType end = min(npoints, source->GetNumberOfPoints());
626         int       dim = min(3,       array->GetNumberOfComponents());
627         for (vtkIdType ptId = 0; ptId < end; ++ptId) {
628           for (int i = 0; i < dim; ++i) p[i] = array->GetComponent(ptId, i);
629           points->SetPoint(ptId, p);
630         }
631         copy = nullptr;
632       } else {
633         if (sourceAttr == pd_as_cd.GetPointer()) {
634           const auto sourceName = (array->GetName()  ? array->GetName()  : string());
635           const auto targetName = (pd[i]._TargetName ? pd[i]._TargetName : string());
636           if (IsCategoricalArrayName(sourceName) || IsCategoricalArrayName(targetName)) {
637             vtkDataArray *labels = GetSourceArray("point", sourcePD, pd[i], case_sensitive);
638             array = ConvertPointToCellLabels(source, labels, label_mode);
639           }
640         }
641         copy = Copy(array, ntuples, pd[i]._TargetDataType);
642       }
643     }
644     if (copy) {
645       if (pd[i]._TargetName) copy->SetName(pd[i]._TargetName);
646       AddArray(targetAttr, copy, pd[i]._TargetAttribute);
647     }
648   }
649   if (point_mask_name) {
650     vtkSmartPointer<vtkDataArray> mask = NewVTKDataArray(VTK_UNSIGNED_CHAR);
651     mask->SetName(point_mask_name);
652     mask->SetNumberOfComponents(1);
653     mask->SetNumberOfTuples(npoints);
654     mask->FillComponent(0, .0);
655     for (vtkIdType ptId = 0; ptId < npoints && ptId < source->GetNumberOfPoints(); ++ptId) {
656       mask->SetComponent(ptId, 0, 1.0);
657     }
658     targetPD->AddArray(mask);
659   }
660 
661   // Add cell data
662   for (size_t i = 0; i < cd.size(); ++i) {
663     if (cd[i]._PointCellConversion) {
664       sourceAttr = cd_as_pd;
665       targetAttr = targetPD;
666       ntuples    = npoints;
667     } else {
668       sourceAttr = sourceCD;
669       targetAttr = targetCD;
670       ntuples    = ncells;
671     }
672     array = GetSourceArray("cell", sourceAttr, cd[i], case_sensitive);
673     if (sourceAttr == cd_as_pd.GetPointer()) {
674       const auto sourceName = (array->GetName()  ? array->GetName()  : string());
675       const auto targetName = (cd[i]._TargetName ? cd[i]._TargetName : string());
676       if (IsCategoricalArrayName(sourceName) || IsCategoricalArrayName(targetName)) {
677         vtkDataArray *labels = GetSourceArray("cell", sourceCD, cd[i], case_sensitive);
678         array = ConvertCellToPointLabels(source, labels, label_mode);
679       }
680     }
681     copy = Copy(array, ntuples, cd[i]._TargetDataType);
682     if (cd[i]._TargetName) copy->SetName(cd[i]._TargetName);
683     AddArray(targetAttr, copy, cd[i]._TargetAttribute);
684   }
685   if (cell_mask_name) {
686     vtkSmartPointer<vtkDataArray> mask = NewVTKDataArray(VTK_UNSIGNED_CHAR);
687     mask->SetName(cell_mask_name);
688     mask->SetNumberOfComponents(1);
689     mask->SetNumberOfTuples(ncells);
690     mask->FillComponent(0, .0);
691     for (vtkIdType cellId = 0; cellId < ncells && cellId < source->GetNumberOfCells(); ++cellId) {
692       mask->SetComponent(cellId, 0, 1.0);
693     }
694     targetCD->AddArray(mask);
695   }
696 
697   // Write resulting data set
698   WritePointSet(output_name, target);
699 }
700