Mercurial > hg > orthanc-stl

comparison Sources/Plugin.cpp @ 29:62abf3c523f9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
big endian support, retrieve axes from referenced volume
author Sebastien Jodogne <s.jodogne@gmail.com>
date Thu, 04 Apr 2024 16:29:01 +0200
parents 410003c50b17
children 3570c23764d4
comparison
equal deleted inserted replaced
28:410003c50b17 29:62abf3c523f9
51 51
52 #include <nifti1_io.h> 52 #include <nifti1_io.h>
53 53
54 #define ORTHANC_PLUGIN_NAME "stl" 54 #define ORTHANC_PLUGIN_NAME "stl"
55 55
56 static const char* const STL_SOP_CLASS_UID = "1.2.840.10008.5.1.4.1.1.104.3";
57
58 56
59 // Forward declaration 57 // Forward declaration
60 void ReadStaticAsset(std::string& target, 58 void ReadStaticAsset(std::string& target,
61 const std::string& path); 59 const std::string& path);
62 60
160 158
161 #include <dcmtk/dcmdata/dcdeftag.h> 159 #include <dcmtk/dcmdata/dcdeftag.h>
162 #include <dcmtk/dcmdata/dcfilefo.h> 160 #include <dcmtk/dcmdata/dcfilefo.h>
163 #include <dcmtk/dcmdata/dcitem.h> 161 #include <dcmtk/dcmdata/dcitem.h>
164 #include <dcmtk/dcmdata/dcsequen.h> 162 #include <dcmtk/dcmdata/dcsequen.h>
163 #include <dcmtk/dcmdata/dcuid.h>
165 164
166 class Extent2D : public boost::noncopyable 165 class Extent2D : public boost::noncopyable
167 { 166 {
168 private: 167 private:
169 bool isEmpty_; 168 bool isEmpty_;
345 double GetZ() const 344 double GetZ() const
346 { 345 {
347 return z_; 346 return z_;
348 } 347 }
349 348
349 double ComputeSquaredNorm() const
350 {
351 return x_ * x_ + y_ * y_ + z_ * z_;
352 }
353
350 double ComputeNorm() const 354 double ComputeNorm() const
351 { 355 {
352 return sqrt(x_ * x_ + y_ * y_ + z_ * z_); 356 return sqrt(ComputeSquaredNorm());
353 } 357 }
354 358
355 void Normalize() 359 void Normalize()
356 { 360 {
357 double norm = ComputeNorm(); 361 double norm = ComputeNorm();
373 377
374 static double DotProduct(const Vector3D& a, 378 static double DotProduct(const Vector3D& a,
375 const Vector3D& b) 379 const Vector3D& b)
376 { 380 {
377 return a.GetX() * b.GetX() + a.GetY() * b.GetY() + a.GetZ() * b.GetZ(); 381 return a.GetX() * b.GetX() + a.GetY() * b.GetY() + a.GetZ() * b.GetZ();
382 }
383
384 static bool AreParallel(const Vector3D& a,
385 const Vector3D& b)
386 {
387 if (IsNear(a.ComputeSquaredNorm(), 1) &&
388 IsNear(b.ComputeSquaredNorm(), 1))
389 {
390 return IsNear(std::abs(Vector3D::DotProduct(a, b)), 1);
391 }
392 else
393 {
394 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange,
395 "Only applicable to normalized vectors");
396 }
378 } 397 }
379 }; 398 };
380 399
381 400
382 401
676 assert(polygons_[i] != NULL); 695 assert(polygons_[i] != NULL);
677 delete polygons_[i]; 696 delete polygons_[i];
678 } 697 }
679 } 698 }
680 699
681 const std::string& GetPatient() const 700 const std::string& GetPatientId() const
682 { 701 {
683 return patientId_; 702 return patientId_;
684 } 703 }
685 704
686 const std::string& GetStudyInstanceUid() const 705 const std::string& GetStudyInstanceUid() const
749 double slicesSpacing_; 768 double slicesSpacing_;
750 double minProjectionAlongNormal_; 769 double minProjectionAlongNormal_;
751 double maxProjectionAlongNormal_; 770 double maxProjectionAlongNormal_;
752 size_t slicesCount_; 771 size_t slicesCount_;
753 772
754 static bool AreNormalsParallel(const Vector3D& a,
755 const Vector3D& b)
756 {
757 return IsNear(std::abs(Vector3D::DotProduct(a, b)), 1);
758 }
759
760 bool LookupProjectionIndex(size_t& index, 773 bool LookupProjectionIndex(size_t& index,
761 double z) const 774 double z) const
762 { 775 {
763 if (slicesCount_ == 0) 776 if (slicesCount_ == 0)
764 { 777 {
815 { 828 {
816 isValid = true; 829 isValid = true;
817 slicesNormal_ = normal; 830 slicesNormal_ = normal;
818 } 831 }
819 832
820 if (AreNormalsParallel(normal, slicesNormal_)) 833 if (Vector3D::AreParallel(normal, slicesNormal_))
821 { 834 {
822 // This is a valid slice (it is parallel to the normal) 835 // This is a valid slice (it is parallel to the normal)
823 const Vector3D& point = structures.GetPolygon(i).GetPoint(0); 836 const Vector3D& point = structures.GetPolygon(i).GetPoint(0);
824 projections.push_back(Vector3D::DotProduct(point, slicesNormal_)); 837 projections.push_back(Vector3D::DotProduct(point, slicesNormal_));
825 } 838 }
1000 bool ProjectAlongNormal(double& z, 1013 bool ProjectAlongNormal(double& z,
1001 const StructurePolygon& polygon) const 1014 const StructurePolygon& polygon) const
1002 { 1015 {
1003 Vector3D normal; 1016 Vector3D normal;
1004 if (polygon.IsCoplanar(normal) && 1017 if (polygon.IsCoplanar(normal) &&
1005 AreNormalsParallel(normal, slicesNormal_)) 1018 Vector3D::AreParallel(normal, slicesNormal_))
1006 { 1019 {
1007 z = Vector3D::DotProduct(polygon.GetPoint(0), slicesNormal_); 1020 z = Vector3D::DotProduct(polygon.GetPoint(0), slicesNormal_);
1008 return true; 1021 return true;
1009 } 1022 }
1010 else 1023 else
1062 } 1075 }
1063 } 1076 }
1064 }; 1077 };
1065 1078
1066 1079
1080 static void WriteFloat(Orthanc::ChunkedBuffer& buffer,
1081 float value,
1082 Orthanc::Endianness endianness)
1083 {
1084 switch (endianness)
1085 {
1086 case Orthanc::Endianness_Little:
1087 buffer.AddChunk(&value, sizeof(float));
1088 break;
1089
1090 case Orthanc::Endianness_Big:
1091 {
1092 uint8_t tmp[4];
1093 tmp[0] = reinterpret_cast<uint8_t*>(&value) [3];
1094 tmp[1] = reinterpret_cast<uint8_t*>(&value) [2];
1095 tmp[2] = reinterpret_cast<uint8_t*>(&value) [1];
1096 tmp[3] = reinterpret_cast<uint8_t*>(&value) [0];
1097 buffer.AddChunk(&tmp, sizeof(float));
1098 break;
1099 }
1100
1101 default:
1102 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
1103 }
1104 }
1105
1106
1107 static void WriteInteger(Orthanc::ChunkedBuffer& buffer,
1108 uint32_t value,
1109 Orthanc::Endianness endianness)
1110 {
1111 switch (endianness)
1112 {
1113 case Orthanc::Endianness_Little:
1114 buffer.AddChunk(&value, sizeof(uint32_t));
1115 break;
1116
1117 case Orthanc::Endianness_Big:
1118 {
1119 uint8_t tmp[4];
1120 tmp[0] = reinterpret_cast<uint8_t*>(&value) [3];
1121 tmp[1] = reinterpret_cast<uint8_t*>(&value) [2];
1122 tmp[2] = reinterpret_cast<uint8_t*>(&value) [1];
1123 tmp[3] = reinterpret_cast<uint8_t*>(&value) [0];
1124 buffer.AddChunk(&tmp, sizeof(uint32_t));
1125 break;
1126 }
1127
1128 default:
1129 throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented);
1130 }
1131 }
1132
1133
1067 static void EncodeSTL(std::string& target /* out */, 1134 static void EncodeSTL(std::string& target /* out */,
1068 vtkPolyData& mesh /* in */) 1135 vtkPolyData& mesh /* in */)
1069 { 1136 {
1070 // TODO - Conversion to little endian on big endian 1137 const Orthanc::Endianness endianness = Orthanc::Toolbox::DetectEndianness();
1071 1138
1072 Orthanc::ChunkedBuffer buffer; 1139 Orthanc::ChunkedBuffer buffer;
1073 1140
1074 uint8_t header[80]; 1141 uint8_t header[80];
1075 memset(header, 0, sizeof(header)); 1142 memset(header, 0, sizeof(header));
1076 buffer.AddChunk(header, sizeof(header)); 1143 buffer.AddChunk(header, sizeof(header)); // This doesn't depend on endianness
1077 1144
1078 uint32_t n = mesh.GetNumberOfCells(); 1145 WriteInteger(buffer, mesh.GetNumberOfCells(), endianness);
1079 buffer.AddChunk(&n, sizeof(n));
1080 1146
1081 for (vtkIdType i = 0; i < mesh.GetNumberOfCells(); i++) 1147 for (vtkIdType i = 0; i < mesh.GetNumberOfCells(); i++)
1082 { 1148 {
1083 vtkCell* cell = mesh.GetCell(i); 1149 vtkCell* cell = mesh.GetCell(i);
1084 vtkTriangle* triangle = dynamic_cast<vtkTriangle*>(cell); 1150 vtkTriangle* triangle = dynamic_cast<vtkTriangle*>(cell);
1091 triangle->GetPoints()->GetPoint(2, p2); 1157 triangle->GetPoints()->GetPoint(2, p2);
1092 1158
1093 double normal[3]; 1159 double normal[3];
1094 vtkTriangle::ComputeNormal(p0, p1, p2, normal); 1160 vtkTriangle::ComputeNormal(p0, p1, p2, normal);
1095 1161
1096 float d[4 * 3] = { 1162 WriteFloat(buffer, normal[0], endianness);
1097 static_cast<float>(normal[0]), static_cast<float>(normal[1]), static_cast<float>(normal[2]), 1163 WriteFloat(buffer, normal[1], endianness);
1098 static_cast<float>(p0[0]), static_cast<float>(p0[1]), static_cast<float>(p0[2]), 1164 WriteFloat(buffer, normal[2], endianness);
1099 static_cast<float>(p1[0]), static_cast<float>(p1[1]), static_cast<float>(p1[2]), 1165 WriteFloat(buffer, p0[0], endianness);
1100 static_cast<float>(p2[0]), static_cast<float>(p2[1]), static_cast<float>(p2[2]) }; 1166 WriteFloat(buffer, p0[1], endianness);
1101 buffer.AddChunk(d, sizeof(d)); 1167 WriteFloat(buffer, p0[2], endianness);
1168 WriteFloat(buffer, p1[0], endianness);
1169 WriteFloat(buffer, p1[1], endianness);
1170 WriteFloat(buffer, p1[2], endianness);
1171 WriteFloat(buffer, p2[0], endianness);
1172 WriteFloat(buffer, p2[1], endianness);
1173 WriteFloat(buffer, p2[2], endianness);
1102 1174
1103 uint16_t a = 0; 1175 uint16_t a = 0;
1104 buffer.AddChunk(&a, sizeof(a)); 1176 buffer.AddChunk(&a, sizeof(a)); // This doesn't depend on endianness
1105 } 1177 }
1106 1178
1107 buffer.Flatten(target); 1179 buffer.Flatten(target);
1108 } 1180 }
1109 1181
1168 1240
1169 return true; 1241 return true;
1170 } 1242 }
1171 1243
1172 1244
1173 bool EncodeStructureSetMesh(std::string& stl, 1245 static Orthanc::ParsedDicomFile* LoadInstance(const std::string& instanceId)
1174 const StructureSet& structureSet, 1246 {
1175 const StructureSetGeometry& geometry, 1247 std::string dicom;
1176 const std::set<std::string>& roiNames, 1248
1177 unsigned int resolution, 1249 if (!OrthancPlugins::RestApiGetString(dicom, "/instances/" + instanceId + "/file", false))
1178 bool smooth) 1250 {
1179 { 1251 throw Orthanc::OrthancException(Orthanc::ErrorCode_UnknownResource);
1180 if (resolution < 1) 1252 }
1253 else
1254 {
1255 return new Orthanc::ParsedDicomFile(dicom);
1256 }
1257 }
1258
1259
1260 static void GetReferencedVolumeAxes(Vector3D& axisX,
1261 Vector3D& axisY,
1262 const StructureSet& structureSet)
1263 {
1264 // This is a rough guess for the X/Y axes
1265 axisX = Vector3D(1, 0, 0);
1266 axisY = Vector3D(0, 1, 0);
1267
1268 // Look for one instance from the referenced volume
1269 const std::string& sopInstanceUid = structureSet.GetPolygon(0).GetReferencedSopInstanceUid();
1270
1271 Json::Value response;
1272 if (OrthancPlugins::RestApiPost(response, "/tools/lookup", sopInstanceUid, false))
1273 {
1274 if (response.type() != Json::arrayValue)
1275 {
1276 throw Orthanc::OrthancException(Orthanc::ErrorCode_NetworkProtocol);
1277 }
1278
1279 bool first = true;
1280
1281 for (Json::Value::ArrayIndex i = 0; i < response.size(); i++)
1282 {
1283 if (response[i].type() != Json::objectValue)
1284 {
1285 throw Orthanc::OrthancException(Orthanc::ErrorCode_NetworkProtocol);
1286 }
1287
1288 if (Orthanc::SerializationToolbox::ReadString(response[i], "Type") == "Instance")
1289 {
1290 if (first)
1291 {
1292 const std::string& instanceId = Orthanc::SerializationToolbox::ReadString(response[i], "ID");
1293 std::unique_ptr<Orthanc::ParsedDicomFile> reference(LoadInstance(instanceId));
1294
1295 std::string imageOrientation;
1296 if (reference->GetTagValue(imageOrientation, Orthanc::DICOM_TAG_IMAGE_ORIENTATION_PATIENT))
1297 {
1298 std::vector<std::string> items;
1299 Orthanc::Toolbox::TokenizeString(items, imageOrientation, '\\');
1300
1301 double x1, x2, x3, y1, y2, y3;
1302
1303 if (items.size() == 6 &&
1304 MyParseDouble(y1, items[0]) &&
1305 MyParseDouble(y2, items[1]) &&
1306 MyParseDouble(y3, items[2]) &&
1307 MyParseDouble(x1, items[3]) &&
1308 MyParseDouble(x2, items[4]) &&
1309 MyParseDouble(x3, items[5]))
1310 {
1311 axisX = Vector3D(x1, x2, x3);
1312 axisY = Vector3D(y1, y2, y3);
1313 }
1314 }
1315 }
1316 else
1317 {
1318 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError,
1319 "Multiple instances with the same SOP Instance UID");
1320 }
1321 }
1322 }
1323 }
1324 }
1325
1326
1327 static bool EncodeStructureSetMesh(std::string& stl,
1328 const StructureSet& structureSet,
1329 const StructureSetGeometry& geometry,
1330 const std::set<std::string>& roiNames,
1331 unsigned int resolution,
1332 bool smooth)
1333 {
1334 if (resolution < 1 ||
1335 structureSet.GetPolygonsCount() < 1)
1181 { 1336 {
1182 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); 1337 throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange);
1183 } 1338 }
1184 1339
1185 if (!IsNear(1, geometry.GetSlicesNormal().ComputeNorm())) 1340 if (!IsNear(1, geometry.GetSlicesNormal().ComputeNorm()))
1186 { 1341 {
1187 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); 1342 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
1188 } 1343 }
1189 1344
1190 // TODO - Axes could be retrieved from the referenced CT volume 1345 Vector3D axisX, axisY;
1191 Vector3D axisX(1, 0, 0); 1346 GetReferencedVolumeAxes(axisX, axisY, structureSet);
1192 Vector3D axisY = Vector3D::CrossProduct(geometry.GetSlicesNormal(), axisX); 1347
1348 Vector3D axisZ = Vector3D::CrossProduct(axisX, axisY);
1193 1349
1194 if (!IsNear(1, axisX.ComputeNorm()) || 1350 if (!IsNear(1, axisX.ComputeNorm()) ||
1195 !IsNear(1, axisY.ComputeNorm())) 1351 !IsNear(1, axisY.ComputeNorm()) ||
1352 !Vector3D::AreParallel(axisZ, geometry.GetSlicesNormal()))
1196 { 1353 {
1197 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); 1354 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError);
1198 } 1355 }
1199 1356
1200 Extent2D extent; 1357 Extent2D extent;
1213 memset(volume->GetScalarPointer(), 0, resolution * resolution * depth); 1370 memset(volume->GetScalarPointer(), 0, resolution * resolution * depth);
1214 1371
1215 for (size_t i = 0; i < structureSet.GetPolygonsCount(); i++) 1372 for (size_t i = 0; i < structureSet.GetPolygonsCount(); i++)
1216 { 1373 {
1217 const StructurePolygon& polygon = structureSet.GetPolygon(i); 1374 const StructurePolygon& polygon = structureSet.GetPolygon(i);
1218 if (roiNames.find(polygon.GetRoiName()) == roiNames.end()) 1375 if (roiNames.find(polygon.GetRoiName()) != roiNames.end())
1219 { 1376 {
1220 // This polygon doesn't correspond to a ROI of interest 1377 // This polygon corresponds to a ROI of interest
1221 continue; 1378
1222 } 1379 size_t z;
1223 1380 if (geometry.LookupSliceIndex(z, polygon))
1224 size_t j; 1381 {
1225 if (!geometry.LookupSliceIndex(j, polygon)) 1382 std::vector<Orthanc::ImageProcessing::ImagePoint> points;
1226 { 1383 points.reserve(polygon.GetPointsCount());
1227 throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); 1384
1228 } 1385 for (size_t j = 0; j < polygon.GetPointsCount(); j++)
1229 1386 {
1230 std::vector<Orthanc::ImageProcessing::ImagePoint> points; 1387 const Vector3D& point = polygon.GetPoint(j);
1231 points.reserve(polygon.GetPointsCount()); 1388 double x = (Vector3D::DotProduct(point, axisX) - extent.GetMinX()) / extent.GetWidth() * static_cast<double>(resolution);
1232 for (size_t j = 0; j < polygon.GetPointsCount(); j++) 1389 double y = (Vector3D::DotProduct(point, axisY) - extent.GetMinY()) / extent.GetHeight() * static_cast<double>(resolution);
1233 { 1390 points.push_back(Orthanc::ImageProcessing::ImagePoint(static_cast<int32_t>(std::floor(x)),
1234 const Vector3D& point = polygon.GetPoint(j); 1391 static_cast<int32_t>(std::floor(y))));
1235 double x = (Vector3D::DotProduct(point, axisX) - extent.GetMinX()) / extent.GetWidth() * static_cast<double>(resolution); 1392 }
1236 double y = (Vector3D::DotProduct(point, axisY) - extent.GetMinY()) / extent.GetHeight() * static_cast<double>(resolution); 1393
1237 points.push_back(Orthanc::ImageProcessing::ImagePoint(static_cast<int32_t>(std::floor(x)), 1394 Orthanc::ImageAccessor slice;
1238 static_cast<int32_t>(std::floor(y)))); 1395 slice.AssignWritable(Orthanc::PixelFormat_Grayscale8, resolution, resolution, resolution /* pitch */,
1239 } 1396 reinterpret_cast<uint8_t*>(volume->GetScalarPointer()) + z * resolution * resolution);
1240 1397
1241 Orthanc::ImageAccessor slice; 1398 XorFiller filler(slice);
1242 slice.AssignWritable(Orthanc::PixelFormat_Grayscale8, resolution, resolution, resolution /* pitch */, 1399 Orthanc::ImageProcessing::FillPolygon(filler, points);
1243 reinterpret_cast<uint8_t*>(volume->GetScalarPointer()) + j * resolution * resolution); 1400 }
1244 1401 }
1245 XorFiller filler(slice); 1402 }
1246 Orthanc::ImageProcessing::FillPolygon(filler, points); 1403
1247 } 1404 volume->SetSpacing(extent.GetWidth() / static_cast<double>(resolution),
1248 1405 extent.GetHeight() / static_cast<double>(resolution),
1249 volume->SetSpacing( 1406 geometry.GetSlicesSpacing());
1250 extent.GetWidth() / static_cast<double>(resolution),
1251 extent.GetHeight() / static_cast<double>(resolution),
1252 (geometry.GetMaxProjectionAlongNormal() - geometry.GetMinProjectionAlongNormal()) / static_cast<double>(depth));
1253 1407
1254 // TODO 1408 // TODO
1255 // volume->SetOrigin() 1409 // volume->SetOrigin()
1256 1410
1257 return EncodeVolume(stl, volume.Get(), resolution, smooth); 1411 return EncodeVolume(stl, volume.Get(), resolution, smooth);
1258 }
1259
1260
1261 static Orthanc::ParsedDicomFile* LoadInstance(const std::string& instanceId)
1262 {
1263 std::string dicom;
1264
1265 if (!OrthancPlugins::RestApiGetString(dicom, "/instances/" + instanceId + "/file", false))
1266 {
1267 throw Orthanc::OrthancException(Orthanc::ErrorCode_UnknownResource);
1268 }
1269 else
1270 {
1271 return new Orthanc::ParsedDicomFile(dicom);
1272 }
1273 } 1412 }
1274 1413
1275 1414
1276 void ListStructures(OrthancPluginRestOutput* output, 1415 void ListStructures(OrthancPluginRestOutput* output,
1277 const char* url, 1416 const char* url,
1510 std::unique_ptr<Orthanc::ParsedDicomFile> dicom(LoadInstance(instanceId)); 1649 std::unique_ptr<Orthanc::ParsedDicomFile> dicom(LoadInstance(instanceId));
1511 DcmDataset& dataset = *dicom->GetDcmtkObject().getDataset(); 1650 DcmDataset& dataset = *dicom->GetDcmtkObject().getDataset();
1512 1651
1513 std::string stl; 1652 std::string stl;
1514 if (GetStringValue(dataset, DCM_MIMETypeOfEncapsulatedDocument) != Orthanc::MIME_STL || 1653 if (GetStringValue(dataset, DCM_MIMETypeOfEncapsulatedDocument) != Orthanc::MIME_STL ||
1515 GetStringValue(dataset, DCM_SOPClassUID) != STL_SOP_CLASS_UID || 1654 GetStringValue(dataset, DCM_SOPClassUID) != UID_EncapsulatedSTLStorage ||
1516 !dicom->GetTagValue(stl, Orthanc::DICOM_TAG_ENCAPSULATED_DOCUMENT)) 1655 !dicom->GetTagValue(stl, Orthanc::DICOM_TAG_ENCAPSULATED_DOCUMENT))
1517 { 1656 {
1518 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadRequest, "DICOM instance not encapsulating a STL model: " + instanceId); 1657 throw Orthanc::OrthancException(Orthanc::ErrorCode_BadRequest, "DICOM instance not encapsulating a STL model: " + instanceId);
1519 } 1658 }
1520 else 1659 else