From dd09f02f53be621f8a3f48e20c6e2ab9563efc98 Mon Sep 17 00:00:00 2001 From: Arfon Smith Date: Fri, 23 Jan 2015 15:25:36 -0600 Subject: [PATCH] Adding some samples to verify new heuristics --- samples/Mathematica/HeyexImport.m | 344 +++++++++++++++++++++++++ samples/Mercury/switch_detection_bug.m | 46 ++++ 2 files changed, 390 insertions(+) create mode 100644 samples/Mathematica/HeyexImport.m create mode 100644 samples/Mercury/switch_detection_bug.m diff --git a/samples/Mathematica/HeyexImport.m b/samples/Mathematica/HeyexImport.m new file mode 100644 index 00000000..457a7e92 --- /dev/null +++ b/samples/Mathematica/HeyexImport.m @@ -0,0 +1,344 @@ +(* Mathematica Package *) +(* Created with IntelliJ IDEA and the Mathematica Language plugin *) + +(* :Title: Importer for the RAW data-format of the Heidelberg Eye Explorer (known as HEYEX) *) + +(* :Context: HeyexImport` *) + +(* :Author: Patrick Scheibe pscheibe@trm.uni-leipzig.de *) + +(* :Package Version: 1.0 *) + +(* :Mathematica Version: 8.0 *) + +(* :Copyright: Patrick Scheibe, 2013-2015 *) + +(* :Discussion: This package registers a new importer which can load the RAW data-format exported by a + Heidelberg Spectralis OCT. The import-functionality can access different information contained + in a file: + 1. The file header which contains meta data like when the patient was scanned etc + 2. The scanned volume data + 3. Images which represent slices of the scanned volume + 4. The Scanning laser ophthalmoscopy (SLO) image which is taken with every scanned patient + 5. The segmentation data for different retina layers provided by the software + +*) + +(* :Keywords: Import, Heyex, OCT, Spectralis, Heidelberg Engineering *) + +BeginPackage[ "HeyexImport`" ] + +HeyexEyePosition::usage = "HeyexEyePosition[file] tries to extract which eye was scanned, left or right."; + +HeyexImport::wrongHdr = "Error importing OCT data. Broken/Wrong file?"; + + +Begin[ "`Private`" ]; + +(* + Registration of all import possibilities for the Heidelberg OCT. +*) + +ImportExport`RegisterImport[ + "Heyex" , + { + "FileHeader" :> importHeader, + { "Data" , n_Integer} :> (importData[n][##]&), + "Data" :> importData, + { "Images" , n_Integer} :> (importImages[n][##]&), + "Images" :> importImages, + "SLOImage" :> importSLOImage, + "SegmentationData" :> importSegmentation, + { "SegmentationData" , n_Integer} :> (importSegmentation[n][##]&), + "DataSize" :> importDataSize, + importData + }, + + { + "Image3D" :> (Image3D["Data" /. #1]&) + }, + + "AvailableElements" -> {"FileHeader", "Data", "DataSize", "Images", "SLOImage", "SegmentationData", "Image3D"} +]; + + +If[Quiet[Check[TrueQ[Compile[{}, 0, CompilationTarget -> "C"][] == 0], False]], + $compileTarget = CompilationTarget -> "C", + $compileTarget = CompilationTarget -> "MVM" +]; + + +(* + Helper function which reads data from a stream. This is + only a unification so I can map the read function over a + list. +*) +read[{id_String, type_String}, str_] := +id -> BinaryRead[str, type]; +read[{type_String, n_Integer}, str_] := BinaryReadList[str, type, n]; +read[{id_String, {type_String, n_Integer}}, str_] := id -> BinaryReadList[str, type, n]; +(* + Note that when reading bytes explicitly I convert them to + a string and remove any zeroes at the end. +*) +read[{id_String, { "Byte" , n_Integer}}, str_] := +id -> StringJoin[ + FromCharacterCode /@ (Rest[ + NestList[BinaryRead[str, "Byte" ] &, Null, + n]] /. {chars___Integer, Longest[0 ...]} :> {chars})]; + +(* + The layout of a file exported with "Raw Export" + + ***************** + * File Header * + ***************** + * SLO Image * + ***************** + * B-Scan #0 * + ***************** + * ..... * + ***************** + * B-Scan #n-1 * + ***************** +*) + +With[{i = "Integer32", f = "Real32", d = "Real64", b = "Byte"}, + + $fileHeaderInfo = Transpose[{ + { + "Version" , "SizeX" , "NumBScans" , "SizeZ" , "ScaleX" , "Distance" , + "ScaleZ" , "SizeXSlo" , "SizeYSlo" , "ScaleXSlo" , "ScaleYSlo" , + "FieldSizeSlo" , "ScanFocus" , "ScanPosition" , "ExamTime" , + "ScanPattern" , "BScanHdrSize" , "ID" , "ReferenceID" , "PID" , + "PatientID" , "Padding" , "DOB" , "VID" , "VisitID" , "VisitDate" , + "Spare" + }, + { + {b, 12}, i, i, i, d, d, d, i, i, d, d, i, d, {b, 4}, {i, 2}, i, i, + {b, 16}, {b, 16}, i, {b, 21}, {b, 3}, d, i, {b, 24}, d, {b, 1840} + } + }]; + + $bScanHeaderInfo = Transpose[{ + { + "Version" , "BScanHdrSize" , "StartX" , "StartY" , "EndX" , "EndY" , + "NumSeg" , "OffSeg" , "Quality" , "Spare" + }, + {{b, 12}, i, d, d, d, d, i, i, f, {b, 196}} + }]; +]; + + +isHeyexRawFormat[{"Version" -> version_String, "SizeX" -> _Integer, "NumBScans" -> _Integer, _Rule..}] /; StringMatchQ[version, "HSF-OCT" ~~__] := True ; +isHeyexRawFormat[___] := False; + +readFileHeader[str_InputStream] := With[{hdr = Quiet[read[#, str]] & /@ $fileHeaderInfo}, + hdr /; TrueQ[isHeyexRawFormat[hdr]] +]; +readFileHeader[___] := (Message[HeyexImport::wrongHdr]; Throw[$Failed]); + + +(* Reads the camera image of the retina. Note that you must have the + information from the fileheader and you must be at the right position + of the file stream for this.*) +readSLOImage[str_InputStream, fileHdr : {(_String -> _) ..}] := + Image[Partition[ + BinaryReadList[str, "Byte" , "SizeXSlo" * "SizeYSlo" /. fileHdr], + "SizeXSlo" /. fileHdr], "Byte" ]; + +skipSLOImage[str_InputStream, fileHdr : {(_String -> _) ..}] := + Skip[str, "Byte" , "SizeXSlo" * "SizeYSlo" /. fileHdr]; + + +(* One single BScan consists itself again of a header and a data part *) +readBScanHeader[str_InputStream, fileHdr : {(_String -> _) ..}] := + Module[{i = "Integer32", f = "Real32", d = "Real64", b = "Byte", + bScanHdr}, + bScanHdr = read[#, str] & /@ Transpose[{ + { "Version" , "BScanHdrSize" , "StartX" , "StartY" , "EndX" , "EndY" , + "NumSeg" , "OffSeg" , "Quality" , "Spare" }, + {{b, 12}, i, d, d, d, d, i, i, f, {b, 196}}} + ]; + AppendTo[bScanHdr, + read[{ "SegArray" , { "Real32" , + "NumSeg" * "SizeX" /. bScanHdr /. fileHdr}}, str] + ]; + (* + This is horrible slow, therefore I just skip the fillbytes + + AppendTo[bScanHdr, + read[{"Fillbytes", {"Byte", + "BScanHdrSize" - 256 - "NumSeg"*"SizeX"*4 /. bScanHdr /. + fileHdr}}, str] + ] + *) + Skip[str, "Byte" , "BScanHdrSize" - 256 - "NumSeg" * "SizeX" * 4 /. bScanHdr /. fileHdr]; + AppendTo[bScanHdr, "FillBytes" -> None] + ] + +skipBScanHeader[str_InputStream, fileHdr : {(_String -> _) ..}] := + Skip[str, "Byte" , "BScanHdrSize" /. fileHdr]; + +readBScanData[str_InputStream, fileHdr : {(_String -> _) ..}] := + Module[{}, + Developer`ToPackedArray[ + Partition[read[{ "Real32" , "SizeX" * "SizeZ" /. fileHdr}, str], + "SizeX" /. fileHdr]] + ]; + +skipBScanData[str_InputStream, fileHdr : {(_String -> _) ..}] := + Skip[str, "Byte" , "SizeX" * "SizeZ" * 4 /. fileHdr]; + +skipBScanBlocks[str_InputStream, fileHdr : {(_String -> _) ..}, n_Integer] := + Skip[str, "Byte" , n * ("BScanHdrSize" + "SizeX" * "SizeZ" * 4) /. fileHdr]; + + +importHeader[filename_String, ___] := Module[ + {str, header}, + str = OpenRead[filename, BinaryFormat -> True]; + header = readFileHeader[str]; + Close[str]; + "FileHeader" -> header +]; + + +(* Imports the dimension of the scanned volume. *) +importDataSize[filename_String, r___] := Module[{header = importHeader[filename]}, + "DataSize" -> ({"NumBScans", "SizeZ", "SizeXSlo"} /. ("FileHeader" /. header)) +] + +importSLOImage[filename_String, ___] := Module[ + {str, header, slo}, + str = OpenRead[filename, BinaryFormat -> True]; + header = readFileHeader[str]; + slo = readSLOImage[str, header]; + Close[str]; + "SLOImage" -> slo +] + +importData[filename_String, ___] := Module[ + {str, header, nx, n, data}, + str = OpenRead[filename, BinaryFormat -> True]; + header = readFileHeader[str]; + {nx, n} = { "SizeX" , "SizeX" * "SizeZ"} /. header; + skipSLOImage[str, header]; + data = Table[ + skipBScanHeader[str, header]; + Partition[read[{ "Real32" , n}, str], nx], + {"NumBScans" /. header} + ]; + Close[str]; + "Data" -> Developer`ToPackedArray[data] +]; + +importData[num_Integer][filename_String, ___] := Module[ + {str, header, nx, n, data}, + str = OpenRead[filename, BinaryFormat -> True]; + header = readFileHeader[str]; + {nx, n} = { "SizeX" , "SizeX" * "SizeZ"} /. header; + skipSLOImage[str, header]; + skipBScanBlocks[str, header, Max[Min["NumBScans" /. header, num - 1], 0] ]; + skipBScanHeader[str, header]; + data = Partition[read[{ "Real32" , n}, str], nx]; + Close[str]; + {"Data" -> {num -> Developer`ToPackedArray[data]}} +]; + +(* + As suggested in the Heidelberg OCT Manual the importer will adjust + the graylevels when importing images. Since this is very time-consuming + for the whole scanned volume, I use an optimized version of this function. +*) +With[{$compileTarget = $compileTarget}, $adjustGraylevelFunc := ($adjustGraylevelFunc = Compile[{{values, _Real, 2}}, + Map[Floor[255.0 * Min[Max[0.0, #], 1.0]^(0.25) + 0.5] &, values, {2}], + RuntimeAttributes -> {Listable}, + Parallelization -> True, + RuntimeOptions -> "Speed", + $compileTarget +])]; + +importImages[filename_String, ___] := Module[ + {data}, + data = "Data" /. importData[filename]; + "Images" -> (Image[#, "Byte" ]& /@ $adjustGraylevelFunc[data]) +] + +importImages[imageNumber_Integer][filename_String, ___] := Module[ + {data}, + data = {imageNumber /. ("Data" /. importData[imageNumber][filename])}; + {"Images" -> {imageNumber -> (Image[#, "Byte" ]& @@ $adjustGraylevelFunc[data])}} +]; + +importSegmentation[filename_String, ___] := Module[ + {str, header, data}, + str = OpenRead[filename, BinaryFormat -> True]; + header = readFileHeader[str]; + skipSLOImage[str, header]; + data = Table[ + Module[{bScanHeader, t}, + {t, bScanHeader} = Timing@readBScanHeader[str, header]; + skipBScanData[str, header]; + bScanHeader + ], {"NumBScans" /. header} + ]; + Close[str]; + (* + The BScanHeaderData contain the segmentation vectors as a single list + of numbers. Before returning the result, I check how many segmentations + there are inside the BScan an I transform the segmentation value list + into separate vectors and call them "ILM", "RPE" and "NFL" like described + in the manual + *) + "SegmentationData" -> Function[{bhdr}, + Block[{numVecs = "NumSeg" /. bhdr, vecNames, nx = "SizeX" /. header}, + If[numVecs > 0, + vecNames = Take[{ "ILM" , "RPE" , "NFL" }, numVecs]; + bhdr /. ("SegArray" -> vec_) :> Sequence @@ (Rule @@@ Transpose[{vecNames, Partition[vec, nx]} ]), + bhdr + ] + ]] /@ data +] + +importSegmentation[num_Integer][filename_String, ___] := Module[ + {str, header, bhdr}, + str = OpenRead[filename, BinaryFormat -> True]; + header = readFileHeader[str]; + skipSLOImage[str, header]; + skipBScanBlocks[str, header, Max[Min["NumBScans" /. header, num - 1], 0] ]; + bhdr = readBScanHeader[str, header]; + Close[str]; + (* See doc above *) + {"SegmentationData" -> {num -> Block[ + {numVecs = "NumSeg" /. bhdr, vecNames, nx = "SizeX" /. header}, + If[ numVecs > 0, + vecNames = Take[{ "ILM" , "RPE" , "NFL" }, numVecs]; + bhdr /. ("SegArray" -> vec_) :> Sequence @@ (Rule @@@ Transpose[{vecNames, Partition[vec, nx]} ]), + bhdr + ] + ] + }} +] + +(* Extracts which eye was scanned. This is stored in the header of the file *) +(* OD stands for oculus dexter which is latin for "right eye" and OS stands + for oculus sinister which is latin for "left eye" *) +HeyexEyePosition[file_String /; FileExistsQ[file]] := Module[{position}, + Check[ + position = "ScanPosition" /. Import[file, { "Heyex" , "FileHeader" }]; + Switch[ + position, + "OD" , + Right, + "OS" , + Left, + _, + $Failed + ], + $Failed + ] +]; + +End[] + +EndPackage[] \ No newline at end of file diff --git a/samples/Mercury/switch_detection_bug.m b/samples/Mercury/switch_detection_bug.m new file mode 100644 index 00000000..f96542c3 --- /dev/null +++ b/samples/Mercury/switch_detection_bug.m @@ -0,0 +1,46 @@ +% This is a regression test for a bug in switch detection +% where it was preferring incomplete switches to complete +% one-case switches, and hence inferring the wrong determinism. + +%------------------------------------------------------------------------------% + +:- module switch_detection_bug. + +:- interface. + +:- type note ---> note(rank, modifier, octave). + +:- type rank ---> c ; d ; e ; f ; g ; a ; b . + +:- type modifier ---> natural ; sharp ; flat . + +:- type octave == int. + +:- type qualifier ---> maj ; min . + +:- pred next_topnote(note, qualifier, note). +:- mode next_topnote(in, in, out) is multi. + +%------------------------------------------------------------------------------% + +:- implementation. + +next_topnote(note(c, _, Oct), _, note(d, natural, Oct)). +next_topnote(note(d, _, Oct), _, note(c, natural, Oct)). +next_topnote(note(d, _, Oct), maj, note(e, natural, Oct)). +next_topnote(note(d, _, Oct), min, note(e, flat, Oct)). +next_topnote(note(e, _, Oct), _, note(d, natural, Oct)). +next_topnote(note(e, _, Oct), _, note(f, natural, Oct)). +next_topnote(note(f, _, Oct), maj, note(e, natural, Oct)). +next_topnote(note(f, _, Oct), min, note(e, flat, Oct)). +next_topnote(note(g, _, Oct), _, note(f, natural, Oct)). +next_topnote(note(g, _, Oct), min, note(a, flat, Oct)). +next_topnote(note(g, _, Oct), maj, note(a, natural, Oct)). +next_topnote(note(a, _, Oct), _, note(g, natural, Oct)). +next_topnote(note(a, _, Oct), min, note(b, flat, Oct)). +next_topnote(note(a, _, Oct), maj, note(b, natural, Oct)). +next_topnote(note(b, _, Oct), maj, note(a, natural, Oct)). +next_topnote(note(b, _, Oct), min, note(a, flat, Oct)). + +%------------------------------------------------------------------------------% +