vtkRAdapter.cxx 20.7 KB
Newer Older
Thomas Otahal's avatar
   
Thomas Otahal committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkRAdapter.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
/*-------------------------------------------------------------------------
  Copyright 2009 Sandia Corporation.
  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
  the U.S. Government retains certain rights in this software.
-------------------------------------------------------------------------*/

#include "vtkObjectFactory.h"
#include "vtkRAdapter.h"
#include "vtkAbstractArray.h"
#include "vtkStdString.h"
#include "vtkStringArray.h"
#include "vtkTable.h"
28
#include "vtkTree.h"
Thomas Otahal's avatar
   
Thomas Otahal committed
29
30
31
32
33
34
35
#include "vtkArray.h"
#include "vtkArrayExtents.h"
#include "vtkArrayCoordinates.h"
#include "vtkTypedArray.h"
#include "vtkVariantArray.h"
#include "vtkDoubleArray.h"
#include "vtkIntArray.h"
36
37
38
#include "vtkDataArrayCollection.h"
#include "vtkArrayData.h"
#include "vtkDataObjectCollection.h"
39
40
41
42
43
44
#include "vtkTreeDFSIterator.h"
#include "vtkSmartPointer.h"
#include "vtkNew.h"
#include "vtkEdgeListIterator.h"
#include "vtkDataSetAttributes.h"
#include "vtkMutableDirectedGraph.h"
Thomas Otahal's avatar
   
Thomas Otahal committed
45

46
#include <map>
Thomas Otahal's avatar
   
Thomas Otahal committed
47
48

#include <stdio.h>
49
#include <cassert>
Thomas Otahal's avatar
   
Thomas Otahal committed
50
51
52
53
54
55
56
57
58

#include "R.h"
#include "Rdefines.h"
#include "R_ext/Parse.h"
#include "R_ext/Rdynload.h"


vtkStandardNewMacro(vtkRAdapter);

59
60
61
namespace
{

Thomas Otahal's avatar
   
Thomas Otahal committed
62
63
64
65
66
67
68
69
70
71
72
int R_FindArrayIndex(vtkArrayCoordinates& coordinates, const vtkArrayExtents& extents)
{

  vtkIdType i;
  int ret = 0;
  vtkIdType divisor = 1;
  vtkIdType d = coordinates.GetDimensions();

  for(i = 0; i < d; ++i)
    {
    ret = ret + coordinates[i]*divisor;
73
    divisor *= extents[i].GetSize();
Thomas Otahal's avatar
   
Thomas Otahal committed
74
75
76
77
78
79
    }

  return(ret);

}

80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
} // End anonymous namespace

//----------------------------------------------------------------------------
vtkRAdapter::vtkRAdapter()
{

  this->vad =  vtkArrayData::New();
  this->vdoc = vtkDataObjectCollection::New();
  this->vdac = vtkDataArrayCollection::New();

}

//----------------------------------------------------------------------------
vtkRAdapter::~vtkRAdapter()
{

  if(this->vad)
    {
    this->vad->Delete();
    }

  if(this->vdoc)
    {
    this->vdoc->Delete();
    }

  if(this->vdac)
    {
    this->vdac->Delete();
    }

}

Thomas Otahal's avatar
   
Thomas Otahal committed
113
114
115
116
117
118
119
vtkDataArray* vtkRAdapter::RToVTKDataArray(SEXP variable)
{

  int i;
  int j;
  int nr;
  int nc;
120
121
122
  vtkDoubleArray * result;
  double * data;

Thomas Otahal's avatar
   
Thomas Otahal committed
123
124
125
126
127
128

  if( Rf_isMatrix(variable) || Rf_isVector(variable) )
    {
    nc = Rf_ncols(variable);
    nr = Rf_nrows(variable);

129
    result = vtkDoubleArray::New();
Thomas Otahal's avatar
   
Thomas Otahal committed
130
131
132
133
134
135
136
137
138
139

    result->SetNumberOfTuples(nr);
    result->SetNumberOfComponents(nc);

    data = new double[nc];

    for(i=0;i<nr;i++)
      {
      for(j=0;j<nc;j++)
        {
140
141
142
143
144
145
146
147
148
149
150
151
        if ( isReal(variable) )
          {
          data[j] = REAL(variable)[j*nr + i];
          }
        else if ( isInteger(variable) )
          {
          data[j] = static_cast<double>(INTEGER(variable)[j*nr + i]);
          }
        else
          {
          vtkErrorMacro(<< "Bad return variable, tried REAL and INTEGER.");
          }
Thomas Otahal's avatar
   
Thomas Otahal committed
152
153
154
155
        result->InsertTuple(i,data);
        }
      }

156
    delete [] data;
157
    this->vdac->AddItem(result);
Thomas Otahal's avatar
   
Thomas Otahal committed
158
    result->Delete();
Thomas Otahal's avatar
   
Thomas Otahal committed
159
160
161
162
163
164
165
166
167
    return(result);
    }
  else
    {
    return(0);
    }

}

168

Thomas Otahal's avatar
   
Thomas Otahal committed
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
SEXP vtkRAdapter::VTKDataArrayToR(vtkDataArray* da)
{

  SEXP a;
  int nr;
  int nc;
  int i;
  int j;
  double* data;

  nr = da->GetNumberOfTuples();
  nc = da->GetNumberOfComponents();

  PROTECT(a = Rf_allocMatrix(REALSXP,nr,nc));

  for(i=0;i<nr;i++)
    {
    for(j=0;j<nc;j++)
      {
      data = da->GetTuple(i);
      REAL(a)[j*nr + i] = data[j];
      }
    }

  return(a);

}

vtkArray* vtkRAdapter::RToVTKArray(SEXP variable)
{

200
201
202
  vtkArray::SizeT i;
  vtkArray::DimensionT j;
  vtkArray::DimensionT ndim;
Thomas Otahal's avatar
   
Thomas Otahal committed
203
204
205
206
207
  SEXP dims;
  vtkArrayExtents extents;
  vtkTypedArray<double>* da;
  da = vtkTypedArray<double>::SafeDownCast(vtkArray::CreateArray(vtkArray::DENSE, VTK_DOUBLE));

208
209
  dims = getAttrib(variable, R_DimSymbol);
  ndim = static_cast<vtkArray::DimensionT>(length(dims));
Thomas Otahal's avatar
   
Thomas Otahal committed
210

Thomas Otahal's avatar
   
Thomas Otahal committed
211
  if (!isMatrix(variable)&&!isArray(variable)&&isVector(variable))
212
213
214
215
    {
    ndim = 1;
    }

Thomas Otahal's avatar
   
Thomas Otahal committed
216
217
  extents.SetDimensions(ndim);

Thomas Otahal's avatar
   
Thomas Otahal committed
218
  if (isMatrix(variable)||isArray(variable))
219
    {
220
    for(j=0;j<ndim;j++)
221
      {
222
      extents[j] = vtkArrayRange(0,INTEGER(dims)[j]);
223
224
225
      }
    }
  else
Thomas Otahal's avatar
   
Thomas Otahal committed
226
    {
227
    extents[0] = vtkArrayRange(0,length(variable));
Thomas Otahal's avatar
   
Thomas Otahal committed
228
229
230
231
232
233
234
235
236
237
238
    }

  da->Resize(extents);

  vtkArrayCoordinates index;

  index.SetDimensions(ndim);

  for(i=0;i<da->GetSize();i++)
    {
    da->GetCoordinatesN(i,index);
239
240
241
242
243
244
245
246
247
248
249
250
    if ( isReal(variable) )
      {
      da->SetVariantValue(index,REAL(variable)[i]);
      }
    else if ( isInteger(variable) )
      {
      da->SetVariantValue(index,static_cast<double>(INTEGER(variable)[i]));
      }
    else
      {
      vtkErrorMacro(<< "Bad return variable, tried REAL and INTEGER.");
      }
Thomas Otahal's avatar
   
Thomas Otahal committed
251
252
    }

253
  this->vad->AddArray(da);
Thomas Otahal's avatar
   
Thomas Otahal committed
254
  da->Delete();
Thomas Otahal's avatar
   
Thomas Otahal committed
255
256
257
258
259
260
261
262
263
  return(da);

}

SEXP vtkRAdapter::VTKArrayToR(vtkArray* da)
{

  SEXP a;
  SEXP dim;
264
265
  vtkArray::SizeT i;
  vtkArray::DimensionT j;
Thomas Otahal's avatar
   
Thomas Otahal committed
266
267
268
269
  vtkArrayCoordinates coords;

  PROTECT(dim = Rf_allocVector(INTSXP, da->GetDimensions()));

270
  assert(da->GetExtents().ZeroBased());
271
  for(j=0;j<da->GetDimensions();j++)
Thomas Otahal's avatar
   
Thomas Otahal committed
272
    {
273
    INTEGER(dim)[j] = da->GetExtents()[j].GetSize();
Thomas Otahal's avatar
   
Thomas Otahal committed
274
275
276
277
278
279
280
281
282
    }

  PROTECT(a = Rf_allocArray(REALSXP, dim));

  for(i=0;i<da->GetSize();i++)
    {
    REAL(a)[i] = 0.0;
    }

283
  assert(da->GetExtents().ZeroBased());
Thomas Otahal's avatar
   
Thomas Otahal committed
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
  for(i=0;i<da->GetNonNullSize();i++)
    {
    da->GetCoordinatesN(i,coords);
    REAL(a)[R_FindArrayIndex(coords,da->GetExtents())] = da->GetVariantValue(coords).ToDouble();
    }

  UNPROTECT(1);

  return(a);

}

SEXP vtkRAdapter::VTKTableToR(vtkTable* table)
{

  SEXP a;
  SEXP b;
  SEXP names;
  int i;
  int j;
  int nr = table->GetNumberOfRows();
  int nc = table->GetNumberOfColumns();
  vtkVariant data;

  PROTECT(a = allocVector(VECSXP, nc));
  PROTECT(names = allocVector(STRSXP, nc));

  for(j=0;j<nc;j++)
    {
    SET_STRING_ELT(names,j,mkChar(table->GetColumn(j)->GetName()));
    if(vtkDataArray::SafeDownCast(table->GetColumn(j)))
      {
      PROTECT(b = allocVector(REALSXP,nr));
      SET_VECTOR_ELT(a,j,b);
      for(i=0;i<nr;i++)
        {
        data = table->GetValue(i,j);
        REAL(b)[i] = data.ToDouble();
        }
      }
    else
      {
      PROTECT(b = allocVector(STRSXP,nr));
      SET_VECTOR_ELT(a,j,b);
      for(i=0;i<nr;i++)
        {
        data = table->GetValue(i,j);
        SET_STRING_ELT(b,i,mkChar(data.ToString().c_str()));
        }
      }
    }

  setAttrib(a,R_NamesSymbol,names);

  return(a);

}

vtkTable* vtkRAdapter::RToVTKTable(SEXP variable)
{

  int i;
  int j;
  int nr;
  int nc;
  SEXP names;
  vtkVariant v;
  vtkTable* result;

  if( isMatrix(variable) )
    {
    nc = Rf_ncols(variable);
    nr = Rf_nrows(variable);
    result = vtkTable::New();
358
359
360
361
362
363
364
365
366
367
368
369
    names = getAttrib(variable, R_DimNamesSymbol);
    if(!isNull(names))
     {
       vtkStringArray* rowNames = vtkStringArray::New();
       for (i = 0; i < nr; ++i)
         {
         std::string rowName = CHAR(STRING_ELT(VECTOR_ELT(names,0),i));
         rowNames->InsertNextValue(rowName);
         }
       result->AddColumn(rowNames);
       rowNames->Delete();
     }
Thomas Otahal's avatar
   
Thomas Otahal committed
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416

    for(j=0;j<nc;j++)
      {
      vtkDoubleArray* da = vtkDoubleArray::New();
      da->SetNumberOfComponents(1);
      if(!isNull(names))
        {
        da->SetName(CHAR(STRING_ELT(VECTOR_ELT(names,1),j)));
        }
      else
        {
        v = j;
        da->SetName(v.ToString().c_str());
        }
      for(i=0;i<nr;i++)
        {
        da->InsertNextValue(REAL(variable)[j*nr + i]);
        }
      result->AddColumn(da);
      da->Delete();
      }
    }
  else if( isNewList(variable) )
    {
    nc = length(variable);
    nr = length(VECTOR_ELT(variable,0));
    for(j=1;j<nc;j++)
      {
      if(isReal(VECTOR_ELT(variable,j)) ||
         isInteger(VECTOR_ELT(variable,j)) ||
         isString(VECTOR_ELT(variable,j)) )
        {
        if(length(VECTOR_ELT(variable,j)) != nr)
          {
          vtkGenericWarningMacro(<<"Cannot convert R data type to vtkTable");
          return(0);
          }
        }
      else
        {
        vtkGenericWarningMacro(<<"Cannot convert R data type to vtkTable");
        return(0);
        }
      }

    result = vtkTable::New();
    names = getAttrib(variable, R_NamesSymbol);
417
418
419
420
421
422
423
424
425
426
427
428
429
    SEXP names2 = getAttrib(variable, R_RowNamesSymbol);
    if (!isNull(names2))
      {
       vtkStringArray* rowNames = vtkStringArray::New();
       for (i = 0; i < nr; ++i)
         {
         std::string rowName = CHAR(STRING_ELT(names2, i));
         rowNames->InsertNextValue(rowName);
         }
       result->AddColumn(rowNames);
       rowNames->Delete();
      }

Thomas Otahal's avatar
   
Thomas Otahal committed
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
    vtkAbstractArray *aa;
    for(j=0;j<nc;j++)
      {
      if(isReal(VECTOR_ELT(variable,j)))
        {
        vtkDoubleArray* da = vtkDoubleArray::New();
        da->SetNumberOfComponents(1);
        aa = da;
        for(i=0;i<nr;i++)
          {
          da->InsertNextValue(REAL(VECTOR_ELT(variable,j))[i]);
          }
        }
      else if(isInteger(VECTOR_ELT(variable,j)))
        {
        vtkIntArray* da = vtkIntArray::New();
        da->SetNumberOfComponents(1);
        aa = da;
        for(i=0;i<nr;i++)
          {
          da->InsertNextValue(INTEGER(VECTOR_ELT(variable,j))[i]);
          }
        }
      else
        {
        vtkStringArray* da = vtkStringArray::New();
        da->SetNumberOfComponents(1);
        aa = da;
        for(i=0;i<nr;i++)
          {
          da->InsertNextValue(CHAR(STRING_ELT(VECTOR_ELT(variable,j),i)));
          }
        }

      if(!isNull(names))
        {
        aa->SetName(CHAR(STRING_ELT(names,j)));
        }
      else
        {
        v = j;
        aa->SetName(v.ToString().c_str());
        }
      result->AddColumn(aa);
      aa->Delete();
      }
    }
  else
    {
    vtkGenericWarningMacro(<<"Cannot convert R data type to vtkTable");
    return(0);
    }

483
  this->vdoc->AddItem(result);
Thomas Otahal's avatar
   
Thomas Otahal committed
484
  result->Delete();
Thomas Otahal's avatar
   
Thomas Otahal committed
485
486
487
488
  return(result);

}

489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
SEXP vtkRAdapter::VTKTreeToR(vtkTree* tree)
{
  SEXP r_tree;
  SEXP names;
  SEXP edge;
  SEXP edge_length;
  SEXP Nnode;
  SEXP node_label;
  SEXP tip_label;
  SEXP classname;


  int nedge, nnode, ntip;

  // R phylo tree is a list of 5 elements
  PROTECT(r_tree = allocVector(VECSXP, 5));
  PROTECT(names = allocVector(STRSXP, 5));

  // traverse the tree to reorder the leaf vertices according to the
  // phylo tree numbering rule;
  // newNodeId is the checkup table that maps a vertexId(starting from 0)
510
  // to its corresponding R tree point id (starting from 1)
511
512
513
514
  vtkIdType leafCount = 0;
  vtkTreeDFSIterator* iter = vtkTreeDFSIterator::New();
  iter->SetTree(tree);
  int nVerts = tree->GetNumberOfVertices();
515
  int *newNodeId = new int[nVerts];//including root vertex 0
516
517
518
  while (iter->HasNext())
    {// find out all the leaf nodes, and number them sequentially
    vtkIdType vertexId = iter->Next();
519
    newNodeId[vertexId] = 0;//initialize
520
521
522
523
524
525
526
527
528
529
    if (tree->IsLeaf(vertexId))
      {
      leafCount++;
      newNodeId[vertexId] = leafCount;
      }
    }

  // second tree traverse to reorder the node vertices
  int nodeId = leafCount;
  iter->Restart();
530
  vtkIdType vertexId;
531
532
533
534
535
536
537
538
539
540
  while (iter->HasNext())
    {
    vertexId = iter->Next();
    if (!tree->IsLeaf(vertexId))
      {
      nodeId++;
      newNodeId[vertexId] = nodeId;
      }
    }

541
  nedge = tree->GetNumberOfEdges();
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
  ntip  = leafCount;
  nnode = nedge - ntip + 1;

  // fill in R variables
  PROTECT(edge = allocMatrix(INTSXP, nedge,2));
  PROTECT(Nnode = allocVector(INTSXP, 1));
  PROTECT(tip_label = allocVector(STRSXP, ntip));
  PROTECT(edge_length = allocVector(REALSXP, nedge));
  PROTECT(node_label = allocVector(STRSXP, nnode));
  INTEGER(Nnode)[0] = nnode;

  int * e = INTEGER(edge);
  double * e_len = REAL(edge_length);

  // fill in e and e_len
  vtkSmartPointer<vtkEdgeListIterator> edgeIterator = vtkSmartPointer<vtkEdgeListIterator>::New();
  tree->GetEdges(edgeIterator);
559
  vtkEdgeType vEdge;
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
  int i = 0;
  vtkDoubleArray * weights = vtkDoubleArray::SafeDownCast((tree->GetEdgeData())->GetArray("weight"));
  while(edgeIterator->HasNext())
    {
    vEdge = edgeIterator->Next();
    e[i]  = newNodeId[vEdge.Source] ;
    e[i + nedge] = newNodeId[vEdge.Target];

    int eNum = tree->GetEdgeId(vEdge.Source, vEdge.Target);
    e_len[i] = weights->GetValue(eNum);
    i++;
    }

  // fill in  Nnode , tip_label and  node_label
  // use GetAbstractArray() instead of GetArray()
  vtkStringArray * labels = vtkStringArray::SafeDownCast((tree->GetVertexData())->GetAbstractArray("node name"));
  iter->Restart();
  while (iter->HasNext())
    {// find out all the leaf nodes, and number them sequentially
    vertexId = iter->Next();
    if (tree->IsLeaf(vertexId))
      {
      vtkStdString lab = labels->GetValue(vertexId);
      SET_STRING_ELT(tip_label, newNodeId[vertexId]-1, mkChar(lab.c_str()));
      }
    else
      {
      vtkStdString lab = labels->GetValue(vertexId);
      SET_STRING_ELT(node_label,newNodeId[vertexId]- ntip - 1, mkChar(lab.c_str())); //the starting id of the internal nodes is (ntip + 1)
      }
    }
591
  iter->Delete();
592

593
  // set all elements
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
  SET_VECTOR_ELT(r_tree, 0, edge);
  SET_VECTOR_ELT(r_tree, 1, Nnode);
  SET_VECTOR_ELT(r_tree, 2, tip_label);
  SET_VECTOR_ELT(r_tree, 3, edge_length);
  SET_VECTOR_ELT(r_tree, 4, node_label);

  SET_STRING_ELT(names, 0, mkChar("edge"));
  SET_STRING_ELT(names, 1, mkChar("Nnode"));
  SET_STRING_ELT(names, 2, mkChar("tip.label"));
  SET_STRING_ELT(names, 3, mkChar("edge.length"));
  SET_STRING_ELT(names, 4, mkChar("node.label"));

  setAttrib(r_tree,R_NamesSymbol,names);

  PROTECT(classname = allocVector(STRSXP, 1));
  SET_STRING_ELT(classname, 0, mkChar("phylo"));
  setAttrib(r_tree, R_ClassSymbol, classname);

612
613
  delete [] newNodeId;

614
615
616
617
618
619
620
621
  UNPROTECT(8);
  return r_tree;
}



vtkTree* vtkRAdapter::RToVTKTree(SEXP variable)
{
622
623
624
625
  int numEdges = 0, numEdgeLengths = 0, numNodes = 0, numNodeLabels = 0;
  int *edge = NULL;
  double *edgeLength = NULL;
  bool haveNodeLabels = false;
626
627
628
629
630
631
632
  vtkTree * tree = vtkTree::New();

  if (isNewList(variable))
    {
    int nELT = length(variable);
    if (nELT < 4)
      {
633
      vtkErrorMacro(<<"RToVTKTree(): R tree list does not contain required four elements!");
634
635
636
      return NULL;
      }

637
638
639
640
641
642
643
644
645
    vtkNew<vtkStringArray> nodeLabels;
    vtkNew<vtkStringArray> tipLabels;
    SEXP listNames = getAttrib(variable, R_NamesSymbol);
    SEXP rTipLabels = NULL;
    SEXP rNodeLabels = NULL;

    // collect data from R.  The elements of the tree list are apparently not
    // guaranteed to be in any specific order.
    for (int i = 0; i < nELT; ++i)
646
      {
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
      const char *name = CHAR(STRING_ELT(listNames, i));

      if (strcmp(name, "edge") == 0)
        {
        SEXP rEdge = VECTOR_ELT(variable, i);
        if (isInteger(rEdge))
          {
          edge = INTEGER(rEdge);
          numEdges = length(rEdge)/2;
          }
        else
          {
          vtkErrorMacro(<<"RToVTKTree(): \"edge\" array is not integer type. ");
          return NULL;
          }
        }

      else if (strcmp(name, "Nnode") == 0)
        {
        SEXP rNnode = VECTOR_ELT(variable, i);
        if (isInteger(rNnode))
          {
          numNodes =  INTEGER(rNnode)[0];
          if (length(rNnode) != 1)
            {
            vtkErrorMacro(<<"RToVTKTree(): Expect a single scalar of \"Nnode\". ");
            return NULL;
            }
          }
        else
          {
          vtkErrorMacro(<<"RToVTKTree(): \"Nnode\" is not integer type. ");
          return NULL;
          }
        }
682

683
684
685
686
      else if (strcmp(name, "tip.label") == 0)
        {
        rTipLabels = VECTOR_ELT(variable, i);
        }
687

688
      else if (strcmp(name, "edge.length") == 0)
689
        {
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
        SEXP rEdgeLength = VECTOR_ELT(variable, i);
        numEdgeLengths = length(rEdgeLength);
        if (isReal(rEdgeLength))
          {
          edgeLength = REAL(rEdgeLength);
          }
        }

      else if (strcmp(name, "node.label") == 0)
        {
        // optional node labels
        rNodeLabels = VECTOR_ELT(variable, i);
        haveNodeLabels = true;
        numNodeLabels = length(rNodeLabels);
        }

      else
        {
        vtkWarningMacro(<<"Unexpected tree element encountered: " << name);
709
710
        }
      }
711
712
713
714

    // Perform some safety checks to make sure that the data extracted from R
    // is sane.
    if (numEdges != numEdgeLengths)
715
      {
716
717
      vtkErrorMacro(<<"RToVTKTree(): "
                    << "edgeLength's size does not match up with numEdges.");
718
719
720
      return NULL;
      }

721
    if (haveNodeLabels && numNodeLabels != numNodes)
722
      {
723
724
725
726
727
728
729
730
731
732
      vtkErrorMacro(<<"RToVTKTree(): node.label's size does not match numNodes.");
      return NULL;
      }

    // Populate node & tip label arrays.  If no such labels were specified by R
    // do so with empty strings.
    nodeLabels->SetNumberOfValues(numNodes);
    if (rNodeLabels != NULL && isString(rNodeLabels))
      {
      for (int i = 0; i < numNodes; i++)
733
        {
734
        nodeLabels->SetValue(i, CHAR(STRING_ELT(rNodeLabels, i)));
735
736
        }
      }
737
    else
738
      {
739
      for (int i = 0; i < numNodes; i++)
740
        {
741
        nodeLabels->SetValue(i, "");
742
743
744
        }
      }

745
746
747
    int numTips = numEdges - numNodes + 1;
    tipLabels->SetNumberOfValues(numTips);
    if (rTipLabels != NULL && isString(rTipLabels))
748
      {
749
      for (int i = 0; i < numTips; i++)
750
        {
751
        tipLabels->SetValue(i, CHAR(STRING_ELT(rTipLabels, i)));
752
753
754
755
        }
      }
    else
      {
756
      for (int i = 0; i < numTips; i++)
757
        {
758
        tipLabels->SetValue(i, "");
759
        }
760
      }
761
762
763
764

    //------------  Build the VTKTree -----------
    vtkNew<vtkMutableDirectedGraph> builder;

765
766
    // Create all of the tree vertices (number of edges +1)
    for(int i = 0; i <= numEdges; i++)
767
768
769
770
      {
      builder->AddVertex();
      }

771
    for(int i = 0; i < numEdges; i++)
772
      {
Zack Galbreath's avatar
Zack Galbreath committed
773
774
      // -1 because R vertices begin with 1, whereas VTK vertices begin with 0.
      vtkIdType source = edge[i] - 1;
775
      vtkIdType target = edge[i+numEdges] - 1;
776
777
778
779
780
781
782
      builder->AddEdge(source, target);
      }

    // Create the edge weight array
    vtkNew<vtkDoubleArray> weights;
    weights->SetNumberOfComponents(1);
    weights->SetName("weight");
783
784
    weights->SetNumberOfValues(numEdges);
    for (int i = 0; i < numEdges; i++)
785
      {
786
      weights->SetValue(i, edgeLength[i]);
787
788
789
790
      }
    builder->GetEdgeData()->AddArray(weights.GetPointer());

    // Create the names array
791
792
793
    // In R tree, the numeric id of the vertice is ordered such that the tips
    // are listed first followed by the internal nodes. Match up this order
    // with the label arrays (tip.label and node.label).
794
795
796
    vtkNew<vtkStringArray> names;
    names->SetNumberOfComponents(1);
    names->SetName("node name");
797
798
    names->SetNumberOfValues(numTips + numNodes);
    for (int i = 0; i < numTips; i++)
799
      {
800
      names->SetValue(i, tipLabels->GetValue(i));
801
      }
802
    for (int i = 0; i < numNodes; i++)
803
      {
804
      names->SetValue(i + numTips, nodeLabels->GetValue(i));
805
806
807
808
809
810
811
812
813
      }
    builder->GetVertexData()->AddArray(names.GetPointer());

    if (!tree->CheckedShallowCopy(builder.GetPointer()))
      {
      vtkErrorMacro(<<"Edges do not create a valid tree.");
      return NULL;
      }

814
815
    // Create the "node weight" array for the Vertices, in order to use
    // vtkTreeLayoutStrategy for visualizing the tree using vtkTreeHeatmapItem
816
817
    vtkNew<vtkDoubleArray> nodeWeights;
    nodeWeights->SetNumberOfTuples(tree->GetNumberOfVertices());
818

819
820
821
822
823
824
825
826
827
828
829
    vtkNew<vtkTreeDFSIterator> treeIterator;
    treeIterator->SetStartVertex(tree->GetRoot());
    treeIterator->SetTree(tree);
    while (treeIterator->HasNext())
      {
      vtkIdType vertex = treeIterator->Next();
      vtkIdType parent = tree->GetParent(vertex);
      double weight = 0.0;
      if (parent >= 0)
        {
        weight = weights->GetValue(tree->GetEdgeId(parent, vertex));
Chris Harris's avatar
Chris Harris committed
830
        weight += nodeWeights->GetValue(parent);
831
832
        }
      nodeWeights->SetValue(vertex, weight);
833
      }
834

835
836
837
    nodeWeights->SetName("node weight");
    tree->GetVertexData()->AddArray(nodeWeights.GetPointer());

838
839
    this->vdoc->AddItem(tree);
    tree->Delete();
840
841
842
843
844
845
846
847
848
849
    return tree;
    }
  else
    {
    vtkErrorMacro(<<"RToVTKTree(): R variable is not a list. ");
    return NULL;
    }
}


Thomas Otahal's avatar
   
Thomas Otahal committed
850
851
852
853
854
void vtkRAdapter::PrintSelf(ostream& os, vtkIndent indent)
{

  this->Superclass::PrintSelf(os,indent);

855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
  if(this->vad)
    {
    this->vad->PrintSelf(os,indent);
    }

  if(this->vdoc)
    {
    this->vdoc->PrintSelf(os,indent);
    }

  if(this->vdac)
    {
    this->vdac->PrintSelf(os,indent);
    }


Thomas Otahal's avatar
   
Thomas Otahal committed
871
872
}